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I. INTRODUCTION 


A. PROBLEM STATEMENT 


The problem of motion stability of ships and other marine vehicles has been 
the subject of extensive studies in the past (Comstock, 1977). Most of these 
studies are with regards the ship’s directional stability in open waters under 
open loop conditions. It is well known that in restricted waters such as canals 
or rivers, although it is possible to have positional stability in principle, in 
reality it is never the case. This is due to the destabilizing effects of the bank 
suction forces and moments. These act always in a destabilizing fashion; i.e., 
after a small initial disturbance they produce a force or a moment which tends 
to increase the ship’s path deviation from nominal (Comstock, 1977). Open 
loop conditions, important as they are, rarely represent real life applications. 
Ships traveling in canals are under closed loop control, typically provided 
through the helmsman or an autopilot. It becomes then important to assess 
the stability of the system under finite disturbances and incorporating some 


modeling of closed loop operations. 


B. OBJECTIVES AND OUTLINE 


This thesis utilizes both linear and nonlinear techniques in order to analyze 


the problem of closed loop dynamic stability of ships in canals. Ship modeling 


is discussed in Chapter II and it follows standard ship maneuvering equations 
(Clayton and Bishop, 1982) augmented by a model for bank suction effects 
(Beck, 1976). A Mariner class ship is selected as a model because of the avail- 
ability of data for its hydrodynamic coefficients (Comstock, 1977), (Bernitsas 
and Kekridis, 1984). A heading control law based on Nomoto’s first-order 
model, coupled with a line of sight guidance law is chosen in order to model 
the fundamental behavior of closed loop conditions. Since, in practice control 
actions are hardly ever instantaneous, appropriate time lags are introduced 
in the feedback. These are modeled using the techniques outlined in (Venne, 
1992). Linear eigenvalue analysis is performed in Chapter III in order to 
reveal parametric stability boundaries. It is established that loss of stability 
occurs when a pair of complex conjugate eigenvalues crosses the imaginary axis 
which results in periodic solutions or limit cycles (Guckenheimer and Holmes, 
1983). A nonlinear analysis based on Hopf bifurcation methods (Chow and 
Mallet-Paret, 1977) and (Hassard and Wan, 1978) is performed in Chapter 
IV. Conclusions derived from this study and some recommendations for fur- 
ther extensions are discussed in Chapter V. The basic programs that were used 


to produce the results are included in the Appendix. 


Il. PROBLEM FORMULATION 


A. SHIP DYNAMICS 


Restricting our attention to the horizontal plane, the mathematical model 
consists of the nonlinear sway (translational motion parallel to the vehicle’s 
longitudinal axis) and yaw (rotational motion about the vertical axis ) equa- 
tions of motion. If we consider a moving Cartesian coordinate frame located 


at the geometric center of the vehicle, Newton’s equations of motion are: 


(1) 


mivu+tur+acgr) = Y, 
Igr +mzeg(0v + ur) =e (2) 


where v and r are the relative sway and yaw velocities of the moving vehicle 
with respect to the water, m is the vehicle’s mass, Jz is its moment of in- 
ertia with respect to the vertical axis, and u is the constant forward speed. 
Since all quantities will be considered as dimensionless in this study, in the 
following we consider u to be equal to one. zg is the longitudinal position 
of the ship’s center of gravity with respect to its centroid, and Y,N represent 
the total excitation sway force and yaw moment respectively. Following stan- 
dard ship maneuvering assumptions, these forces can be expressed as the sum 
of quadratic drag terms and first-order memoryless polynomials in v and r, 


which represent added mass and damping due to water. In this way they can 


be expressed by: 


1 1 
Ve You ae 5 Yow” year 


Z 
] 
5 Yeo” +Y,r+Y(v,y) + Y66 , (3) 
I 1 
N = Nyv+Ner+ Nyv t+ 5 Newt + 5 Nerrur® aie 
] 
5 Nerv" +N,r + N(v,y) + Noé , (4) 


where Y,, Nz represent the partial derivatives with respect to the indicated 
variable a and 6 is the effective rudder angle. Y(w,y), N(w,y) represent the 
interaction/proximity forces and moments that arise due to the presence of 


the canal, and shallow water effects. 


The resulting nonlinear differential equations can be non-dimensionalized 
with respect to the constant forward speed of the ship, u and its length, 1. 
The dimensionless time variable is then equal to tu/l. A standard inertial sys- 
tem (x,y) is introduced here where the z-axis points in the assumed nominal 
straight line path and the y-axis is the distance from this nominal path, see 
Figure 1. We assume that the nominal straight line path is along the center- 
line of the canal. In cases where the concept of a geometric centerline is not 
applicable, we assume that the nominal path is along the zero bank suction 
location. ‘This is consistent with recommended navigation practices that are 


currently in use. 


The complete model is presented by the following equations: 


i Hi 
(m — Yj)v — (Y% — mzg)r = Yu + g Yow” i vekvT ee 
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Figure 1: Vehicle geometry and definitions of symbols 


Yr +(¥,-—m)r +Y(b,y) + Yo6 , (5) 
—(N; —mzg)v + (I, - N;)t = Nyv+ = Now ae = Nurrvr?+ 

5 rut’ +(N, — mag)r +N(b,y) + No6 , (6) 
ae (7) 
y=siny +vcosy , (8) 


where the expressions for the ship’s yaw rate as well as its inertial position rate 
have been added, and w is the local heading angle. The interaction/proximity 


forces and moments are expanded to include up to third order terms, 


1 ] 
Y (vy, y) a Yyw a Yyy — a Yow ce 5 Yow” , (9) 
1 . I 3 


D 


as explained in more detail in the folowing section. 


B. HYDRODYNAMIC COEFFICIENTS 


We chose a Mariner class ship as a representative model. Its hydrodynamic 
coefficients and geometric properties were taken from (Comstock, 1977) and 
(Bernitsas and Kekridis, 1984). Results from (Beck, 1976) were used in order 
to model the bank suction forces and moments. Typical results are shown in 
Figure 2. These show force and moment coefficients versus ship deviation (7), 
and for canal width w = 0.4L. Force and moment coefficients are nondimen- 
sionalized with respect to the water density and the ship’s speed and length, 
as is customary in ship maneuvering. Since the suction force and moment 
must change their sign as either y or ~ changes its sign, they must be modeled 
by odd power polynomials, as was done in the previous section. Numerical 
values for the coefficients Y,, Yyyy, Ny, and Nyyy can be found by curve fitting 
of the results of Figure 2. Using a depth to draft ratio of 1.9 we were able to 


estimate, 


Ye e— ee 
Yyyy = 0.468 , 

Ni = 0 00708 
Nyy = 0 


The value of Ny was estimated by “discretizing” the ship in two segments, 


4 
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FIGURE 13 
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Figure 2: Forces and moments due to canal [Beck (1976)] 
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at the bow and the stern, and calculating the suction forces on each part 


individually. Their resultant moment produced, 
Ny = 0.01. 


The value of Nyy was taken to be zero, due to lack of reliable data. The value 
of Yy was estimated to be equal to —Y, = 0.014 which is true for motions along 
the centerline of the canal (Comstock, 1977). Again, due to lack of reliable 


data we took Yyy, = 0. 


C. CONTROL LAW 


The linearized set of equations (5) through (7) can be expressed in the 


following form: 


D5 
Vv = ayy t+ aigu + aygr + aigy + 5}6 , (11) 
rT = agp + aggu + ao3r + agay + 526 , 


where the coefficients a;;, b; are functions of the vehicle hydrodynamic coefh- 


cients and geometric properties and are given below: 














ae = (dg — No ee eee nae 
a2 = ame — N;)Yy + (¥% -— mzg)N,| , 
ag = [Uz - Ni) - m) + (¥%s~ mag)(N,— mae)], 
ne 5 dp = NAY, HOG =areeiNele 


UT 


1 














an = 5—[(Ne~maa)¥p + (m—¥s)NvI, 
ao9 = = [((Na — marg)Yy + (m — Y5)Ny] , 
eae alin Caer an YN. — mee) . 
a aaln —mzg)Yy+(m—Y5)Ny] , 
a 5rltz — N;)¥5 + (Y¥; — mxg)Nz] , 
ae sale Oy Ae 


Dy, = (m—Y%5y)Iz— Ne) — (Ns — mzg)(Y¥; — mz@) . 


A control law that could model a human operator should not be based on 
feedback of the side slip velocity v. Instead, it is more likely that human 
operators will respond to changes in the ship’s heading angle w and rate of 
change of heading, r. Therefore, we choose to base our control law on Nomoto’s 


model, which follows by assuming negligible sway velocity v, 


r=ar+ep+bé, ar) 
where the coefficients are 
a aed Q23 ? 
b= bo ’ 
C = a9) . 


A linear control law can be introduced in the form, 


bo = ki(h — We) + kar , (13) 


where w, 1s the commanded heading angle and the gains kj, ke are computed 
such that the closed loop system (7), (12), and (13) has the desired dynamics. 
The existence of the difference (w — w,) in the control law (13) has the effect of 
trying to point the ship’s longitudinal axis towards the desired heading. The 


characteristic equation of the system is obtained from (7), (12), and (13) as 
s* — (a+ bko)s — (c+ bki) = 0. 
If the desired characteristic equation is 
S- ee  ee 


the control gains can be computed from 


ky = oe 
b b 
os aa 


The natural frequency w, and damping ratio ¢ are selected based on general 
properties of second order systems. Relatively high values of w, and low values 
of ¢ will correspond to a responsive operator, while the opposite is true for 


combinations of low w, and high ¢. 


Finally, in order to take into account the effect of rudder saturation, the 


commanded rudder angle is given by 





6 
§ = Seat tanh ( : | (14) 


sat 


where ég is the slope of 6 at the origin given by (13), and dgat is the saturation 
limit on 6, typically around 0.4 radians. The hyperbolic tangent function is 


used instead of a hard saturation function, because of its differentiability. 
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D. GUIDANCE SCHEME 


Since the previous control law stabilizes the ship to any commanded head- 
ing angle, it must be coupled with an appropriate orientation guidance scheme 
to provide path keeping along the desired track. ‘The simplest guidance scheme 
which models some fundamental aspects of helmsman behavior is a pure pur- 
suit guidance where the commanded heading angle equals the line of sight 
angle, 


te = —tan7? (=) . (15) 


Ld 
as shown in Figure 1. The ship is located at (z,y) and attempts to point its 
longitudinal axis towards a target point which is located ahead of the ship 
on the reference path at a constant preview distance zg. Pursuit guidance is 
achieved by commanding a heading angle w, equal to the line of sight angle 


(15). 


The positional error information y in equation (15), is assumed to lag the 


actual error y by an amount of T seconds, in other words, 
y=yt—-T). (16) 


In this equation, the time lag T models the necessary time that it takes for 
the helmsman to process his path deviation and initiate appropriate corrective 


actions. 


1] 





Ill. LINEAR ANALYSIS 


In this Chapter, we present a linearized analysis of the equations of motion. 
The purpose of this analysis is to assess stability or instability of the equations 
in response to small deviations from the straight line reference track. No 
attempt will be made here to characterize the transient response of the system. 


This can be accomplished through numerical simulations. 


A. COMBINED SYSTEM 


If we incorporate the interaction/proximity forces (9) and (10) in the ship 
dynamics model, equations (5) through (8), we end up with the combined 
system, 

i=in, 

a 3, | 2,1 2 

(m — Y,)v — (Y¥; — mag)r = Yyv + Glow + 5 oer UT + Biro cs 

iF 1 
RY, - m)r Set ele Uden wie? =e Pas + Y56, 


6 6 
1 1 
—(Nz a MZG)v 25 Ge cae N;)r = Nyv a 5 Newr® AP 5 Norevr + (17) 
1 ] 
=Neyyru? + (N, — mag)r + Ny + Nyy + a Newt? + 


Naw" + N66 , 
jy=siny +ucosy. 


Study of the asymptotic properties of this system is the subject of this and 


the following chapters. 
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B. LINEARIZATION 


The linearized form of equations (17) is the following, 
par, 
(m — Y;,)0 =e (Y; =r mzg)r = Te at (be er m)r+ 
Yyw + ae. + ¥56 , 
—(N; — mzg)v + Uz, — N;)f = Nyv + (N, — mzg)r+ 
Nyw + Nyy + N66 , 
y= 


The rudder angle 6 has one of the forms that are developed below. The time 


(18) 


delay operator can be expressed in terms of its Taylor series expansion, 


1 il 2 
y(t-T)=y—Ty+ ae = A eee, (19) 


Practical computations can be performed by truncating equation (19) to first, 


second, or third order. 


1. First order approximation in y 








We have, 

Yt ee 
Or 

yt—-T)=y-T(yty). (20) 
In its linear form, 
6 6 

tanh ( | = - ; 

and, therefore, 
0 = One 


Furthermore, in its linear form equation (15) can be written as 


OH aarti) 
Id 


Y= 


14 


If we incorporate (20) into (13), we derive the linearized first order approxi- 
mation in 6, 
kT 


k 
6 = ky + kor + —y —- ——(v +0). (21) 
Da Lad 


2. Second order approximation in y 


Keeping the second order terms in y(t — T) we get, 
ee 
yt—-T)=y—Ty+ oT iy. (22) 


If we incorporate (22) into (13) along with the linear equations 6 = 69 and 


We = —y(t —T)/xa, we get 











k Ta ie 
6 = kyh + kor + Sy — ——(b + v) + ——(r +0). (23) 
Zd Ld Zan 
3. Third order approximation in y 
In this case, 
ee ea) 
eater SF oy = ty (24) 
and 
ky kit iis ie 
Oy a ee es oe) ot (r+0)- (r+). (25) 
oe on 2a 624 


4. First order approximation in 6 


If a time lag, T, exists on 6 instead of simply y, then 


6 
é — Oc tanh ( : : 
Osat 





where, 


6; = 6o(t —-T) = 59 — Tb. 
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This equation models a time lag associated with the entire application of the 
necessary corrective control action and not just its positional error. Using the 
above equations along with the linear equations 6 = 6; and 7%, = —y/zqg, we 
get 


ky kyT : 
=kypb+kor + ae —k,Tr- ae + v) —koTr . (26) 
d d 


5. First order approximation in both y and 6 
In a more general setting, we can assume that a time lag, 7), exists in the 
control law 6, and a different time lag, T>, is present in the processing of the 


positional error y. Assuming a first order approximation for both, we have 
61 = 6o(t — Ti) = 50 — Tido 


and 
y(t — To) = y — doy . 


Therefore, in this case using the above equations along with the linear equa- 








tions 6 = 6; and w, = —y(t — To)/xq ,we obtain 
k ae 
5 = yy + kor + y — ki —(b +0) - 
kyT: Br 
—=(b +0) + _ 2 eae ee (27) 
d 


Using one of the above expressions, the linearized equations of motion can 


be written as 


b a 

b= ayy + ayqu + aygr + ayay + 16 , (28) 
Tr = aqW + aggu + ao3r + aoay + bod , 

y=ornu, 


where all coefficients have been previously defined. 
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C. SERIES EXPANSIONS OF TIME LAG 


1. First order approximation in y 


In this case we have, 
y(t-T)=y-Ty. 


Substitution of (21) into (28) yields the following matrix system, 








ib OO aro 7 
Cae eee le 22 eee ets v (29) 
r Agi A321 A331 A341 eal 
y 1 1 0 0 y 
where, 
b,k1T 
Ajai = 2 ——— 
Ld 
bk, T 
A221 = a 
oa 
A231 = 413 + bike 
| bik 
A241 = ayy 
a 
bok T 
oie Ce 0k E 
Ld 
bo kyT 
A321 = ag2- ; 
Ld 
A331 = 4a93 + bok 
bok 
A341 = Ce 
Ld 


2. Second order approximation in y 


In this case we have, 


1 
So) ie Ty 


7 


Substitution of (23) into (28) yields the following matrix system, 











i 0 oOo 1 0 db 
v | _ | Aai2 A222 A232 Aga,2 Uv (30) 
r A312 A322 A332 A342 T lies 
y li 1 0 0 y 
where, 
G1) b1k12Xq4 a 6,k,T 
Al,2  — => == 
Ld — 0.96,k iT 
Ghee 
A? — ae 
Id — 0.96, k,T 
yi = Q13Z%_d =r bi korg a 0.50, kT? 
ee, Ld — 0.50,k1T? 
A — y4Zq + by ky 
me pa — 0.5b1 kT? 
bokiT = bok T? 
Azi2 = G21 + bok; — aes ————— A) > 
Ld 20 R 
bokiT bok, T? 
A322 = a22- ee A222 
Ga 2a 
bokiT? bok, T? 
A332 = a3 + boke + ae ae A23,2 
Xd ie 
bok, bok T? 
A342 = ee A242 
Ld 22-4 


3. Third order approximation in y 


A third order approximation in y is based on, 


] 1 a 
Ua) ere ier eae 


Similar algebraic steps as in the previous two approximations result in the 


following eigenvalue problem, 


10 0 0 0 w) 
et i 0 
0 0 B333 0 B353 p |= 
0 0" SO ae y 
0 0 Bs33 0 Bs5,3 v2 
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0 0 
0 0 


1 0 0 
0 0 1 


A3i3 A323 A333 A343 A35,3 


1 i 


0 0 0 


Asi3 As23 As33 A543 A553 


2 3 ~S 


v2 


(31) 


where, v1 = v, vg is the rate of change of v, and the entries of the generalized 


eigenvalue problem (31) are given bellow. Higher order approximations in T 


can not produce any usable stability results, since the B matrix in (31) becomes 


singular. 


A313 
A323 
A33,3 
A343 
A353 
A513 
A592,3 
A533 
A543 
A553 
B33,3 
B35.3 


Bs3.3 














b,kiT 
Qj] + bik, = aaeeeee 
Lg 
b4k1T 
aj2— 
Id 
bi kT? 
Q13 + bike + 
Zia 
bik 
ae Ji*l 
Ld 
bykyT? ; 
20 
bokiT 
Q9] + bok, = 
Id 
bok1T 
a 
Id 
bok iT? 
ao3 + boke + 
Das 
bok 
Toyee wed 
Ld 
bok iT? 
One 
bok T? 
624 
B33.3 
bok1T? 
ie aa 
62g 


bok\T? 
624 





Bs53 = 


4. First order approximation in 6 








Substitution of (26) into (28) yields the following matrix system, 
co OT o- 0 at ab 
QO 1 Boz, 0 v | _ | Ana A224 A234 Ara v (32) 
0 0 Bs34 0 r A314 A324 A334 Asa Tle 
oo 0 ee 1 1 0 0 y 
where, 
b,k\T 
Agia == Gy) 0[ hi" ae 
Ld 
b1k4T 
Aga = 412 — a 
La 
Ag34 = a13 + by kq — b1kiT 
bik 
Ao4 = a4 + —— 
La 
bokyT 
A314 = a1 + bgky — 
Ld 
bokyT 
Azjzo4 = a22 — pees 
Iq 
A334 = a3 + boke — bokiT 
bok 
Azga4 = Goat a 
Lq 
Bo34 = ObykeT 
B33.4 a 55 bokoT 
5. First order approximation in both y and 6 
Substitution of (27) into (28) yields the following matrix system, 
Ler a0 0 Of a 0 0 1 0 yp 
O Boos Bo3z5 0 v | _ | Aas A225 A235 Aras v (33) 
0 Bso5 B335 0 r Agis A325 Ag35 Aga al 
eyo e i iis ia) | Geo aa y 
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where, 


b4k1T 6, kT: 
Aas Cee — i i 
aaa a 
bykiT, = ky T2 
A 5 a2 - —— - -——— 
Ld Ld 
bk, 7477: 
A23.5 €13 + bike — byk1T) + ——— 
A 
bik 
A245 oo —— 
Ld 
bokiT bokiT: 
Asis teeooks = ela oe 
Ld Ld 
b5k1T, bokiT> 
A32.5 2 - + 
“Ld “Ld 
boky TT: 
A335 ao3 + bok — bokiT, + eo 
A 
bok 
A345 ig te 
Ld. 
61k, 71 T: 
ae eee 
Ld 
B35 by koTy 
bok1 TT: 
Peay eee 2 
Ld 
aes 1 + bokoT 


D. EXACT COMPUTATION OF TIME LAG 


The previous analysis using Taylor series expansions of y(t — T’) breaks 
down for approximations beyond third order, as the corresponding generalized 
eigenvalue problem becomes singular. Thus, in order to obtain an exact com- 
putation of the stability curves and check the validity of the calculations, a 


different technique will be performed in this section. This technique is based 


zal 


on frequency response methods and it utilizes Nyquist’s criterion for stability 
[Friedland (1986)].We write the system of equations (11) along with equations 


y=vrtvandd=kwt+koer t+ ae — T), in the Laplace domain, where 
y(t a i) ae ye"® ) (34) 
is the time delay operator. The characteristic equation of the system is, 


1 
8 +015 +959 asst oes ee Bo)—e"* =0, (35) 
d 


where 
QQ, = 412 — 493 — bgke , 
a2 = —@ 14 + 12093 — ao9b) ke — a29a13 + a12bek2 — boky — aay , 
Q3 = —G24 — 013024 — A24b, ko + ao3014 + a14bok — bik 1a02 — 
411422 + bok1a12 + 21412 , 
Q4 = 91014 + 12094 — 422014 — A11G24 + oki 014 — b1k1024 , 
Bg = —bik,, 
By = —a@13b9k, + by k1a93 — bok, , 
Bo = ayqboky — by kya99q — bgkya11 + bi ki a9 , 


The characteristic equation is written as, 


I 
1+—G(s)=0, 
Ld 


where 


(Bos* + Bis + Boje 7 


Cl) =" eee 
(s) s4t+aj,s?+aos* +a3s +04 


(36) 


oe 


is the open loop transfer function , and = denotes the effective position gain. 
For stability, we can utilize the Nyquist criterion which states that the critical 


value of zg can be computed from, 
] . 
—KEGw)I =p for arg G(jw) = —7 (3% 
d 


where 7 denotes the imaginary unit. The argument, arg G(jw), of(36) is, 


' Ss Byw 24 ga aw? 
G =—-wl +t ne —_ 38 
arg G(jw) Ww an By — Bow? Serer es, (38) 

The set of equations (37) result in the critical visibility distance, 

De 2\2 
> as 
ta = ee (Bo — Bow*) (39) 
(a3w — ayw?)? + (w4 — aqw? + a4)? 

The value of w in(38) such that argG(jw) = —7 is the value of the phase 


crossover frequency. After solving for the phase crossover frequency, the critical 
value of xg is obtained by direct substitution of this value of w into equation 


(39). 


E. RESULTS AND DISCUSSION 


Typical results are presented in Figures 3 through 6. All results shown 
are nondimensional unless otherwise stated. The critical value of xg versus wy, 
for zero time lag, and parametrized for different values of the damping ratio 
¢ is shown in Figure 3. Stability is ensured for values of xg greater than its 
critical value. It can be seen that lower values of w, require higher values 


of zg for stability. This means that a less responsive helmsman will need a 
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2 25 S 3.0 4 4.5 5 
omega n 


Figure 3: Critical zg versus wy, for Tg = 0 and various values of ¢ 





2 25 3 3.9 4 4.5 5 
omega n 


Figure 4: Critical zg versus wy for ¢ = 0.8 and various values of T9 
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Figure 5: Critical xg versus wy, for ¢ = 0.8, ‘> = 5sec and various values of 1) 


:| —— With Channel:Effects 
:| —— Without Channel Effects: 
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omega n 


Figure 6: Critical xg versus wy for ¢ = 0.8 and zero time lag: Channel effects 





longer lookahead distance for a stable operation. Similar conclusions hold for 
variations in the damping ratio ¢. In this case, higher values of ¢ correspond 


to better helmsman response which requires less lookahead distance. 


The effect of time lag T2 is shown in Figure 4. Time lags are in seconds. It 
can be seen that reasonable amounts of time lag do not have a serious effect on 
stability, at least in a linear sense. Of course, an amount of time lag may have 
a serious effect on the transient response of the system as well as its ability to 
reject an external disturbance. This can only be established by a systematic 
series of numerical simulations. Similar conclusions hold for non-zero time 


lags T,, as shown in Figure 5. 


Finally, the severe destabilizing effect of the canal is demonstrated by the 
results of Figure 6 where a comparison between canal and open water results is 
presented. It can be seen that an order of magnitude increase in the lookahead 
distance may be required if the same control parameters are to be used in both 


open water and a canal. 


26 





IV. NONLINEAR ANALYSIS 


A. INTRODUCTION 


It can be numerically confirmed that in all cases of stability loss of the pre- 
vious chapter, one pair of complex conjugate eigenvalues of the corresponding 
eigenvalue problem crosses transversally the imaginary axis. A situation like 
this in which a certain parameter is varied such that the real part of one pair 
of complex conjugate eigenvalues of the linearized system matrix crosses zero, 
results in the system leaving its steady state in an oscillatory manner. This 
loss of stability is called Hopf bifurcation and generically occurs in one of two 
ways, supercritical or subcritical. In the supercritical case, stable limit cycles 
are generated after the nominal straight line motion loses its stability. The 
amplitudes of these limit cycles are continuously increasing as the parameter 
distance from its critical value is increased. For small values of this criticality 
distance the resulting limit cycle is of small amplitude and differs little from 
the initial nominal state. In the subcritical case, however, stable limit cycles 
are generated before the nominal state loses its stability. Therefore, depend- 
ing on the initial conditions it is possible to diverge away from the nominal 
straight line path and converge towards a limit cycle even before the nominal 
motion loses its stability. This means that in the subcritical Hopf bifurcation 
case the domain of attraction of the nominal state is decreasing and in fact 


it shrinks to zero as the critical point is approached. Random external dis- 


at 


turbances of sufficient magnitude can throw the vehicle off to an oscillatory 
steady state even though the nominal state may still remain stable. After the 
nominal state becomes unstable, a discontinuous increase in the magnitude of 
motions is observed as there exist no simple stable nearby attractors for the 
vehicle to converge to. Distinction between these two qualitatively different 
types of bifurcation is, therefore, essential in design. The computational pro- 
cedure requires higher order approximations in the equations of motion and is 


the subject of the next section. 


B. DETAILED CALCULATIONS 


The nonlinear equations of motion are written as, 


pom rc, (40) 
0 = ay tayutaygr t+ ayy t 16 , (41) 
F = agp + agqu + ag3r + aogy + bod , (42) 
y = Siny+vcosy~, (43) 


where the contro] law is, 


t—T 
eta tan? tite) (44) 
Ld 


and, including saturation, 





6 
5° = 6.04 tanh ( ; . (45) 
Osat 
We perform a Taylor series expansion of the equations, keeping the first non— 


vanishing nonlinear terms. Due to port/starboard symmetry in the problem 
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this means expansion to third order terms, 


siny =~ — iv, cosy =1— $y°, (46) 
] 
§ =6 -—>-8, 
30%, a0) 
6=kyw+kor + Lie —T)- A134 — T) (48) 
Ld oa 


where for consistency, the term y(t — T) is approximated by its first order 


expansion in 7’, 
y(t-T)=y-Ty=y—-Te+eTv’ —Tv+3Tvy’. (49) 


Substitution of equations (46) to (49) into equations (40) to (45), produces 
the system, 


x = Ax +(x) , (50) 
where the state variables vector is, _ 
> SS na yl] 3 


A is the linearized matrix at equilibrium, and g(x) contains all third order 


terms. 


If T is the matrix of eigenvectors of A evaluated at the critical point zg, 


the linear change of coordinates, 
x=Tz, z2=T'x, (40) 
transforms system (50) into its normal coordinate form, 


2=TUATz+T 'g(Tz) . (52) 


Zo 


At the Hopf bifurcation point, matrix T~!AT takes the form, 


0 —Wpo CAG 

=) a W9 0 0 0 
EE ene ae 

0 0 ~ [Oaea 


where wo is the imaginary part of the critical pair of eigenvalues, and the 
remaining two eigenvalues p and q are negative(or complex conjugate with 
negative real parts) . Therefore in the new system of coordinates (51) the 
dynamics of (50) are governed by a reduced two-dimensional system z, and 
zg since the coordinates z3 and z4 correspond to the eigenvalues p and g and 
are asymptotically stable. For values of zg close to the bifurcation point zg, 


matrix T~'AT is, 


a’e —(wo + w’e) 0 0 

ie Ol etioeeere ave 0 0 

ee 0 0 pipe 0 : 
0 0 0 qtqé 


where, € denotes the criticality difference, 


= (53) 


critical ? 


w = derivative of the real part of the critical eigenvalues with respect to « , 
a = derivative of the imaginary part of the critical eigenvalues 

with respect toe, 
p = derivative of p with respect toe , 


q = derivative of g with respect toe. 


Due to continuity, the eigenvalues p+ p’e and g+q’e remain negative for small 


nonzero values of «. Therefore, the coordinates z3, zq4 correspond to negative 


30 


eigenvalues and are asymptotically stable. Center manifold theory predicts 
that the relationship between the critical coordinates z,, zg and the stable 
coordinates z3, 24 is at least of quadratic order. In fact, due to the symmetry 


in our problem the relationship is cubic, 
21 = O(23,24), 22 = O(23, 24) . 


It follows that because of this order, z3, z4 do not influence the third order 
Taylor series expansions, and can be dropped from the equations. Therefore, 
we can write the reduced system.that describes the essential dynamics of (52) 
in the form, 

Zz) = WEz — (wo tw'e)zo+ Fy(z, 22) , (54) 


Zq = (wo twe)zy + ez + Fo(z}, 20) , (55) 
where F\, Fy contain the third order terms, 


3 2 2 3 
Tider 2422 Ie z2o I T1425 (56) 


Fy (2, 22) 


5 2 2 3 
Fo(21, z2) 7912) + 17222122 + 1932125 +1242 - (57) 


The coefficients r;; are computable from the previous Taylor expansions, and 


they are given below. 


The nonlinear expansion coefficients r;; that are utilized in the definition 


of the cubic stability coefficient K are given by the following: 


"11 3221 + 213£3) + 214l41 + 212d11 + 1213421 , 


T12 119009 + 123339 + Nial42 + N12d12 + 213422 , 
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13 


14 


ill 


22 


123 


24 


where we denote, 


For numerical purposes the critical eigenvector of T must be normalized so 


that its first nonvanishing coefficient is identically equal to 1. The coefficients 


€,; are given by, 


far 
by 


2 2 2 2 
+ OappyMj{ M41 os Obst a a dyur™51M31 =e byrr™M21M3, 


nyi2lo3 + 1313033 + n14l43 + 212413 + 213d93 
ny2o4 + 21334 + Ny4lag + N129d14 + 113d04 
noola + 2331 + noglar + Needy + no3do1 
nozloe + nozlzo + noalge + N22d12 + n23d20 
n22lo3 + n2333 + Nealgs + Noed13 + n23d93 


nooloa + 19334 + No4laa + Noodi4 + no3doq 


Serle Te || eee 


} 


} 


} 


; 
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; 


= byypyMy Mar + SywvM11Mg, + Syyr™M7,M31 + OyrrMyM3, 


Z 2 2 2 
+OyvyM51M41 aI OvyyM21™M 41 te OrryM3,™Ma1 ae OryyM31™M 4, 


+ Our 11M21™M31 + OpyyM11M21M41 + OyryM11M41M31 


3 3 3 3 
+Oyry™M21™41M31 -— OpyyMy, “ir OyyyId) = OrrrM 51 air Oyyyl4y : 


+ byyr(m?,m32 ca 2m 44m12M3}) a Swrr(my2m3, ale 2m 11™31M32) 
+ Sywy(M7M49 + 2m 11m 42™M41) + Sapyy(M4,M12 + 2m 4m41M42) 


2 2 
+d yyr(™51M32 + 2m31M21M22) + burr(mM22M3, + 2m31Mm32M21) 
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= = S yyy (MIM + 2m11M12M21) + Syoo(mM1gM3, + 217m 12M29) 


log 


boa 


+ byyy(m5, Mag + 2m41mM21mM2g2) + Suyy(moom3, + 2maymagmon) 

+ Spry(m3 M42 + 2m41M31™M32) + Spyy(m32m3, + 2mM41:m42M31) 

+ byur[™11™21M32 + m31(m11M22 + M12M21)] 

+ byuy|M11M21M 42 + M41(M11M22 + M12M21)| 

+ 6yrylM11M41M32 + M31(M11M42 + M12™M4))| 

+ byry[M21M41M3q + ™31(mMo1mM429 + M22M41)| 

+d yyy (3miymiz) ala byyy(3m3,m29) ay brrr(3m3,™32) a byyy(3m4ym42) 
Syyu(Migmar + 2miymM12M22) + Syoo(M11M}q + 2m1gMI1M22) 
+5yyr(migmai + 2mi1mM12M39) + Syrr(M11M3q + 2m12mM31M39) 

+ bypy(Minma + 2m44m12M42) + Syyyy(Migmy + 2mj2mM41M42) 
+5yur(mM3omsg1 + 2m21mM29M32) + Surr(Mo1M}q + 2m31m22M39) 
+Syyy(meomay + 2m21m22mM42) + Svyy(Mam4o + 2m41M42M22) 

+ Spry (mogmay + 2m42™m31M32) + Spyy(M31M4q + 2mM41mM42M32) 

+ dyurlm12M22™M31 + M32(m11M22 + M12M21)| 

+byuy|M12M22Ma41 + M42(M11M22 + M12M21))] 

+byrylm12M42M31 + m32(mM11Ma42 + M12M41)| 

+6yry[m22m4gmM31 + M32(Ma1Ma42 + M22M41)] 

+ dyyy(3miim{y) ar byvv(3m21M5p) ih brrr(3m31M 3p) as byyy(3m41m 4p) 
Supe {gM22 + SyvyM1gM39 + SyyrMgM32 + OyrrM12M 3p 

+ SypyM2om a2 + SyyyM1gQM 49 + SyyrM3om 32 + SyrrM22M 5p 


Zz. 2 2 2 
+ Oyyy™M9M 42 ar Ovyy™M22M 49 ar OrryM39M 42 air OryyIN32M 49 
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+6 pyr 127227232 + Oyyy7N12M22M42 + Oyry™M12M42M32 + SyryM22M42M32 


+b yyy te omnia “tn OES ar b yyy 40 , 


fa) en 
bg 
so bn 
boy 
(33, 228 
bg 
sa Ea 
8 
a I 
fai = 7 (max = smu] 
1 1 
49 = -mMi11 (ramos ate 531122 ss sigma: | 
f43 = —M12 (mimes ae smn} a =mizmit | 
_ 1 
l44 = ee ae Ge + sm) u 


The coefficients d;; are given by, 


























dy = = fers — ere err 
d2 = = leraiz — Ne) Pr cone = mag) 
aj3 = = lei32 — Ne) eg ee) 
dig = = [cia(Zz — Nz) + coa(¥e — mag)]|, 
dj) = = [e11(Ns — mag) + cai(m — Yz)] , 
dag = = [ci2(Ng — mag) + c22(m — Y5)] , 
do3 = = [ei3(Vz — mag) + co3(m — Yo)] , 
das = F—Leua(Ne- mo) + exa(m - Yi) 
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where 


C11 


C12 


C13 


C14 


C21 


C22 


C23 


C24 


1 1 1 
gYowm, ar 5 Yurrmg Mar ar 9 -ayils7M31 


I ee 3 
6 pyy~myy + 6 yyy 4l 


iE 1 
2 2 
= vyyv9M 51M 22 = urr (3,9 a 2™m31M32M9)) 


6 2 

1, 9 1 2 

5 Yrev(mam32 + 2m31mM21M22) + 5 yy3M71M12 
1 

5 vvy3M 4,749 

1 il 

5 Yuwdmagmo1 ain 5 Yorr(m3gmor si 2M31M32M22) 
1 1 

5 Yrvu(mogmsi + 2m32M21M22) + 6 pb BMin™M1] 
G say ST gon , 

1 


ee 2 l 2 
g veo ae 9 err 139722 a 5 troviNg9M32 


1 3 1 3 
gr eeem™i2 aif 67 yy 42 » 


] 3 ica. 9 1 9 

g Vow Mar al perc oye ale 5 NrevMg1M31 

I 3 I 3 

aN eww  + gu a1 | 

1 Z I 2 

g Newdmg1 M22 ai: 5 Nore (m31mM22 =o 2313241 ) 
1 2 1 2 

5 Nrwwlmaimse an 2m31M21M22) ala 6 yp3™My1™M12 
1 

5 Ny aMap 

1 2 1 Z 

GV pops ggMa1 = 5 Nurr(m3gmor =a 2313222) 
1 9 iF 2 

5 Nrww(magms + 2m3omM21M22) + g vowed gmail 
1 

= Nuwy3™ag™a1 

1 3 ih 9 1 9 

GV 22 ai 5 NurrM 32M 29 ag 5 Nrovggm 32 

i 3 1 o 

aN gdwMyg + ZNyyy™ go - 
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The coefficients in a third order expansion of the control law are defined 
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C. NONLINEAR COEFFICIENT K 
If we introduce polar coordinates in the form, 
Ze COsiee co — SiN ge (58) 


we can use (54) and (55) to produce an equation describing the rate of change 


of the radial coordinate R, 
R=o'eR+P(O)R? , (59) 
and a similar equation in the rate of change of the angular coordinate @, 
6 = wy +w'e + O(6)R? . (60) 


The system of equations (59) and (60) contains two variables, one of which, 


R, is slowly varying in time, whereas the other one, @ is a fast variable. Then, 


3f 


equation (59) can be averaged over one cycle in 6 to produce an equation with 


constant coefficients and similar stability properties, 
R=o'eR+KR*, (61) 


where, 


1 27 
K= = | P (0) dO = & (3r11 +113 + 722 + 32a) . (62) 


Nontrivial equilibrium solutions of (61) correspond to limit cycles in the 


original system. The nontrivial equilibrium of (61), Ro, is given by, 


Ro =f. (63) 


In our problem, since the trivial equilibrium gains its stability for rg > Ta....01; 
the coefficient a’ is always negative. Therefore, it can be seen from (63) that 
a limit cycle will exist provided K and « have the same sign. The stability 


properties of this limit cycle can be determined by linearization of (61) around 


Ro. The linearized system is, 


R= -—2a’e(R— Ro) , (64) 


and its eigenvalue is, 


B = —2a’e. (65) 


We can see that the Floquet exponent (65) is negative if ¢ is negative, which 
means that K must be negative. Therefore, location and stability of lmit 
cycles depends entirely on the cubic coefficient K which is computable from 


(62). We can summarize our findings as, 
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Figure 7: Nonlinear coefficient K versus wn for dsat = 0.4 and various values of ¢ 


e If K < 0 then limit cycles exist for e < 0 (tq < 2@ ) and they are 


critical 


stable. 


e If K > 0 then limit cycles exist for ¢ > 0 (tq > 2g and they are 


critical ) 


unstable. 


D. RESULTS AND DISCUSSION 


Typical results in terms of the nonlinear stability coefficient K are shown in 
Figures 7 through 10. The general conclusion from this series of graphs is that 
the bifurcations are all strongly subcritical. This means that any linearized 


stability results should be viewed with extreme caution. Limit cycles exist 
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Figure 8: Nonlinear coefficient K versus wy, for ¢ = 0.8 and various values of dsat 
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Figure 9: Nonlinear coefficient K versus wy, for ¢ = 0.8, dsazp = 0.4, T1 = 0, and 
various values of T> (sec) 
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Figure 10: Nonlinear coefficient K versus wp, with and without channel effects 


before stability in the linear sense is lost and a self sustained oscillation in the 
system may develop as a result of an external disturbance even if the nominal 
equilibrium state is still stable. The bifurcations become more subcritical: 
1.e., K is more positive, for smaller values of ¢, as Figure 7 shows. Figure 8 
demonstrates the effect that the saturation limit 6, has on the value of K. 
Higher values of 6,4, although are not related to the critical value of zg result 
in significant changes in the value of K. The general trend is that higher 6a, 
results in less subcritical bifurcations as evidenced by the overall decrease of 
kK. Figure 9 shows the effect of non-zero time lag on the nonlinear stability 
coefficient. It can be seen that, like the linear stability results, the effect is 


minimal. Finally, Figure 10 shows the canal effect on K. This figure was 


Al 





produced for zero time lag, ¢ = 0.8, and ésa: = 0.4. It can be seen that the 
existence of a canal causes the bifurcations to be much more subcritical than 
the open water case. This demonstrates the severe destabilizing effect that the 


canal introduces in both the linearized and the nonlinear analysis. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


This thesis presented a comprehensive study of linear and nonlinear stability 
properties of straightline motion of surface ships in confined waters. The 
classical maneuvering equations of motion incorporating canal suction effects, 
were coupled with appropriate navigation, gauidance, and control laws in order 
to mimic the helmsman’s behavior. The main conclusions from this study can 


be summarized as follows: 


1. There exists a critical preview distance for straightline positional stabil- 


ity. For values less than the critical distance, the system is unstable. 


2. Including canal effects, the critical preview distance may be an order of 
magnitude higher than in open water. If the same preview distance is to 
be used in both cases, the corresponding control law for canal maneuver- 


ing must be considerably more responsive than in open waters. 


3. The critical preview distance is monotonically decreasing for increasing 
control law responsiveness. This means that in order to accomodate 
smaller values of the preview distance, more responsive control laws are 


required . 


4. Physically realizable time lags do not seem to have a significant effect on 


the values of the preview distance. 


5. As the preview distance becomes less than its critical value, one pair of 
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complex conjugate eigenvalues of the linearized system matrix crosses 
the imaginary axis. This corresponds to a bifurcation to periodic so- 
lutions (Hopf bifurcation) and the system exhibits oscillatory behavior, 


also known as limit cycles. 


6. Higher order approximations in the equations of motion were utilized in 
order to assess the stability of the resulting limit cycles. It was found 
that in all cases, the limit cycles were unstable. This has the following 


implications: 


(a) It is possible for the system to lose its stability even before the critical 
preview distance is crossed. This means that all linearized stability 


analysis results should be viewed with extreme caution. 


(b) As the critical preview distance is crossed, it is expected that the 
system will develop limit cycles of large amplitude. This clearly 
presents a dangerous situation which should be avoided in practice, 


by appropriate changes in the design parameters. 


Recommendations for further research include the following: 


1. Simulation studies in order to verify the limit cycle behavior. 


2. Study of different ship characteristics and canal geometry. 
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APPENDIX 


The following is a list and description of the computer programs used in this 
thesis. The programs are written in FORTRAN. Complete printouts of the 
programs follow. Standard eigenvalue/eigenvector numerical analysis subrou- 


tines are required for all programs. 


e THESIS1.FOR 


Calculation of critical zg. First order approximation for time lag T». 


e THESIS2.FOR 


Calculation of critical zz. Second order approximation for time lag T». 


e THESIS3.FOR 


Calculation of critical zz. Third order approximation for time lag T». 


e THESIS4.FOR 


Calculation of critical xg. First order approximation for time lag 7}. 


e THESISS.FOR 
Calculation of critical xg. First order approximation for time lags T; and 


T9. 


e HOPF.FOR 
Calculation of the nonlinear cubic stability coefficient K. Requires the 


output of one of the previous programs as its input. 
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PROGRAM THESIS1.FOR (Time Delay-1st Order Approx. T2) 
BIFURCATION ANALYSIS 


PROGRAM THESIS1 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION K1,K2,L,NR,NV,NDELTA,NPSI,NY,IZ,MASS, 


& NRDOT , NVDOT 


V 


& 


DIMENSION A(4,4) ,FV1(4) ,1V1(4) ,ZZZ(4,4) ,WR(4) ,WI(4) 


OPEN (11,FILE=’BIF1.RES’ ) 
OPEN (12,FILE=’BIF2.RES’) 
OPEN (13,FILE=’BIF3.RES’ ) 
open (15,file=’eig.res1’) 


ehicle Parameters: 
Zz =0.0 

L =528 

RHO =1.94 

XG =0.0 

MASS =0.0088 

U =1.0 


YRDOT = 0.00000 
YVDOT =-0.00912 


YR =+0.00456 
YV =-0.01434 
YPSI = 0.01400 
so 4 = 0.02000 


YDELTA= 0.00278 
NRDOT =-0.00115 
NVDOT = 0.00000 


NR =-0.00296 
NV =-0.00460 
NPSI = 0.01000 
NY =-0.00250 


NDELTA=-0.00139 

WRITE (*,1001) 

READ (*,*) WNMIN,WNMAX, IWN 
WRITE (*,1002) 

READ (*,*) XDMIN , XDMAX , IXD 
WRITE (*,1003) 

READ (*,*) ZETA 

WRITE(*, 1100) 

READ (*,*) TL 

TL=TL*U/L 


DVR =(IZ-NRDOT) * (MASS-YVDOT) - 
(MASS*XG-YRDOT) * (MASS *XG-NVDOT ) 
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& 


& 


AA11=((1IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 

AA12=((IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV)/DVR 

AA21=( (NVDOT-MASS*XG) *YPSI+(MASS-YVDOT) *NPSI) /DVR 

AA22=( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 

AA13=((IZ-NRDOT) * (YR-MASS) + 
(YRDOT~MASS*XG) * (NR-MASS*XG) ) /DVR 

AA23=( (NVDOT-MASS*XG) * (YR-MASS) + 
(MASS-YVDOT) * (NR-MASS*XG) ) /DVR 

AA14=((IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 

AA24=( (NVDOT-MASS*XG) *YY+(MASS-YVDOT) *NY) /DVR 

BBi =((IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 

BB2 =((MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 


ANOM=AA23 
BNOM=BB2 
CNOM=AA21 
EPS =1.D-5 
ILMAX=1500 


DO 1 I=1,IWN 
WRITE (*,2001) I,IWN 
WN=WNMIN+ (I-1) * (WNMAX-WNMIN) / (IWN-1) 
Ki=- (WN**2.D0O/BNOM) ~ (CNOM/BNOM) 
K2=- (ANOM+2.D0*ZETA*WN) /BNOM 
Boe? Jai, TxD 
XD=XDMIN+ (J-1) * (XDMAX-XDMIN) / (IXD-1) 


CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23, 
AA24,BB1,BB2,A,TL) 

CALL RG(4,4,A,WR,WI,0,ZZZ,1IV1i,FV1, IERR) 
CALL DSTABL(DEOS,WR,WI ,FREQ) 

write (15,*)DEOS,XD,WN 


PET Cl-Gh.1)) GO 10 10 
DEOSOQO=DEOS 

XDOO =XD 

LL=0 

COMrGr 2 

DEOSNN=DEOS 

XDNN =XD 
PR=DEOQSNN*DEOSOO 

IF (PR.GT.O.DO) GO TO 3 
LL=LL+1 

TP e<LE.Gi.3) STGP 1000 
IL=0 

XDO=xD00 

XDN=XDNN 

DEOSO=DEOSOO0 


A7 


DEOSN=DEOSNN 

6 XDL=XDO 
XDR=XDN 
DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR) /2.DO 


CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21, 

& AA22,AA23,AA24,BB1,BB2,A,TL) 
CALL RG(4,4,A,WR,WI,0,ZZZ,1V1,FV1, IERR) 
CALL DSTABL(DEOS,WR,WI,FREQ) 
DEOSM=DEOS 
XDM=XD 
PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 
IF (PRL.GT.0.DO) GO TO 5 
XDO=XDL 
XDN=XDM 
DEOSO=DEOSL 
DEOSN=DEOSM 
IL=IL+1 
IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDL-XDM) 
IF (DIF.CGT.EPS) GO TO 6 
XD=XDM 
GO TO 4 

5 IF (PRReGE20.D0) STOP 3200 
XDO=XDM 
XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
IL=IL+1 
IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDM-XDR) 
IF (DIF2GT EPS) GO TO 6 


XD=XDM 
t LLL=10+LL 

WRITE (LLL,*) XD,WN 
3 XDOO=XDNN 


DEOSOO=DEOSNN 
2 CONTINUE 
1 CONTINUE 


1001 FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF WN (dimensionless) ’) 
1002 FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF XD (dimensionless) ’) 
1003 FORMAT (’ ENTER DAMPING RATIO’) 

1100 FORMAT (’ ENTER TIME LAG TL (dimensional) ’) 

2001 FORMAT (215) 


48 


'?) 


?) 


END 
SUBROUTINE DSTABL(DEOS, WR, WI ,OMEGA) 
EVALUATES THE EIGENVALUE WITH THE MAXIMUM REAL PART 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4) ,WI(4) 
DEOS=-1.0D+20 
DO 1 I=1,4 
IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR (TI) 
IJ=1 
CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 


SUBROUTINE LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14, AA21, 


& AA22,AA23,AA24,BB1,BB2,A,TL) 


FORMS THE LINEARIZED MATRIX A (time delay ist order approximation 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION K1,K2 
DIMENSION A(4,4) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 
A(2,1)=AA11+BB1*K1-BB1i*K1*TL/XD 
A(2,2)=AA12-BB1*K1*TL/XD 

A( 2,3) =AA13+BB1*K2 

A(2,4) =AA14+BB1*K1/XD . 
A(3,1)=AA21+BB2*K1-BB2*K1*TL/XD 
A(3,2)=AA22-BB2*K1*TL/XD 
A(3,3)=AA23+BB2*K2 
A(3,4)=AA24+BB2*K1/XD 
A(4,1)=1.0D0 

A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

RETURN 

END 


AQ 


Q 
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PROGRAM THESIS2.FOR (Time Delay-2nd Order Approx T2) 
BIFURCATION ANALYSIS 


PROGRAM THESIS2 

IMPLICIT DOUBLE PRECISION (A-B 0-2) 

DOUBLE PRECISION K1,K2,L,NR,NV,NDELTA,NPSI,NY,IZ,MASS, 
& NRDOT , NVDOT 

DIMENSION A(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4,4) , WR(4) , WI (4) 


OPEN (11,FILE=’BIF1.RES’ ) 
OPEN (12,FILE=’BIF2.RES’) 
OPEN (13,FILE=’BIF3.RES’) 
open (15,file=’eig.resi’) 


Vehicle Parameters: 


TZ =0.0 
L =528 
RHO £=1.94 
XG =0.0 
MASS =0.0088 
U =1.0 


YRDOT = 0.00000 
YVDOT =-0.00912 


XR =+0.00456 
YV =-0.01434 
YPSI = 0.01400 
rey = 0.02000 
YDELTA= 0.00278 


NRDOT =-0.00115 
NVDOT = 0.00000 
NR =-0 .00296 
NV =-0 .00460 


NPSI = 0.01000 
NY =-0 .00250 


NDELTA=-0 .00139 
WRITE (*, 1001) 


READ (*,*) WNMIN, WNMAX, IWN 
WRITE (*,1002) 
READ (*,*) XDMIN, XDMAX, IXD 
WRITE (*,1003) 


00 


READ (*, *) ZETA 
WRITE (*, 1100) 

READ (*,*) TL 
TL=TL*U/L 


DVR =(IZ-NRDOT) * (MASS-YVDOT) - 
& (MASS *XG-YRDOT) * (MASS *XG-NVDOT) 


AA11=((IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12=( (IZ-NRDOT) * YV-(MASS*XG-YRDOT) *NV) /DVR 
AA21=( (NVDOT-MASS*XG) *YPSI+(MASS-YVDOT) *NPSI)/DVR 
AA22=( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13=((IZ-NRDOT) * (YR-MASS) + 


& (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23=((NVDOT-MASS*XG) * (YR-MASS) + 
& (MASS-YVDOT ) * (NR-MASS*XG) ) /DVR 


AA14=( (IZ-NRDOT) * YY+(YRDOT-MASS*XG) *NY) /DVR 
AA24=((NVDOT-MASS*XG) * YY+ (MASS-YVDOT) *NY) /DVR 


BB1 =((IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2 =((MASS-YVDOT) *NDELTA- (MASS *XG-NVDOT) *YDELTA) /DVR 


ANOM=AA23 
BNOM=BB2 
CNOM=AA21 


EPS =1.D-5 
ILMAX=1500 


DO 1 I=1,IWN 
WRITE (*,2001) I,IWN 
WN=WNMIN+(I-1) * (WNMAX-WNMIN) / CIWN-1) 
Ki=- (WN**2.D0/BNOM) -(CNOM/BNOM) 


K2=- (ANOM+2 .DO*ZETA*WN) /BNOM 


DO 2 J=1,IXD 
XD=XDMIN+(J-1) * (XDMAX-XDMIN) / (IXD-1) 


CALL LINEAR(K1 ,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23, 
& AA24 ,BB1,BB2,A,TL) 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL(DEOS ,WR,WI,FREQ) 

write(15,*)DEOS,XD,WN 


ieee Gih. 1) G0. 10. 10 


DEOSOO=DE0S 
XDOO =XD 


ol 


10 


LL=0 

GO TO 2 

DEOSNN=DEOS 

XDNN =XD 
PR=DEOSNN*DEOSOO 

IF (ekeGr.0.D0) GO 10-3 
LL=LL+1 

IF (LE.GT.3) Stor 1000 
ae 

XDO=XD00 

XDN=XDNN 

DEOSO=DE0S00 
DEOSN=DEOSNN 

XDL=XDO 

XDR=XDN 

DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR) /2.D0 


CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21, 
AA22,AA23,AA24,BB1,BB2,A,TL) 

CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1,IERR) 

CALL DSTABL(DEOS,WR,WI, FREQ) 

DEOSM=DE0S 

XDM=XD 

PRL=DEOSL*DEOSM 

PRR=DEOSR*DEOSM 

IF (PRL.GT.0.DO) GO TO 5 

XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 

DIF=DABS (XDL-XDM) 

IF (DIF.GT.EPS) GO TO 6 

XD=XDM 

GO TO 4 

IF (PRR.GT.0.D0) STOP 3200 

XDO=XDM 

XDN=XDR 

DEOSO=DEOSM 

DEOSN=DEOSR 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDM-XDR) 

IF (DIF.GT.EPS) GO TO 6 
XD=XDM 


O2 


QQ 


Q 


— 
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1002 
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1100 
2001 


& 


LLL=10+LL 
WRITE (LLL,*) XD,WN 
XDOO=XDNN 
DEOSOO=DEOSNN 
CONTINUE 
CONTINUE 


FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF WN (dimensionless) ’) 
FORMAT (’? ENTER MIN, MAX, AND INCREMENTS OF XD (dimensionless) ’) 
FORMAT (’ ENTER DAMPING RATIO’) 

FORMAT (’ ENTER TIME LAG TL (dimensional) ’) 

FORMAT (215) 

END 


SUBROUTINE DSTABL(DEOS , WR, WI , OMEGA) 
EVALUATES THE EIGENVALUE WITH THE MAXIMUM REAL PART 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4) ,WI(4) 
DEOS=-1 .0D+20 
DO 1 I=1,4 
IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR (I) 
EL 
CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 


SUBROUTINE LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21, 
AA22 ,AA23,AA24,BB1,BB2,A,TL) 


FORMS THE LINEARIZED MATRIX A (time delay ist order approximation) 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DOUBLE PRECISION ki,k2 

DIMENSION A(4,4) 

RCIG1)=0. ODO 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0DO 
A(2,1)=(AA11*XD+BB1*K1*XD-BB1*K1*TL) / (XD-0 .50D0*BB1*K1*TL*TL) 
A(2,2)=(AA12*XD-BB1*K1*TL) / (XD-0.50D0*BB1*K1*TL*TL) 
A(2,3)=(AA13*XD+BB1*K2*XD+0. 50D0*BB1*K1*TL*TL) / 


& (XD-0.50D0*BB1*K1*TL*TL) 


A(2,4)=(AA14*XD+BB1*K1) / (XD-0 .50D0*BBi*K1*TL*TL) 


o3 


A(3,1)=AA21+BB2*K1-BB2*K1*TL/XD+ (BB2*K1*TL*TL/ (2.0D0O*XD) )*A(2, 1) 
A(3,2)=AA22-BB2*K1*TL/XD+ (BB2*K1*TL*TL/ (2.0D0*XD) ) *A(2, 2) 
A(3,3) =AA23+BB2*K2+ (BB2*K1*TL*TL/ (2.O0D0*XD) )+ 

& (BB2*K1*TL*TL/ (2.0*XD))* A(2,3) 
A(3,4)=AA24+BB2*K1/XD+(BB2*K1*TL*TL/ (2.0D0*XD) ) *A(2,4) 
A(4,1)=1.0D0 

A(4,2)=1.0D0 

A(4,3)=0.0DO 

A(4,4)=0.0DO 

RETURN 

END 


o4 


qa 


PROGRAM THESIS3.FOR (Time Delay-3rd Order Approx T2) 
BIFURCATION ANALYSIS 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION Ki,K2,L,NR,NV,IZ,MASS ,NDELTA,NPSI ,NY, 


& NRDOT , NVDOT 
DIMENS@OieA (59) 5,545,959) ,FVi(5) ,TV1 (59 222(5, 5) ALFR(5), 
& ALPICS) ,BETACS) .WRCS)oWI(5) 


OPEN (11,FILE=’BIF1.RES’) 
OPEN (12, FILE=*BIF2.RES’ ) 
OPEN (13,FILE=’BIF3.RES’ ) 


Vehicle Parameters: 


IZ =0.0 

L =528 
RHO =1.94 
XG =0.0 
MASS =0.0088 
U =1.0 


YRDOT = 0.00000 


MYDOT =——0 00912 
YR =+0 .00456 
XV =-0.01434 
YPSI = 0.01400 
YY = 0.02000 
YDELTA= 0.00278 
NRDOT =-0.00115 
NVDOT = 0.00000 
NR =-0.00296 
NV =-0.00460 
NPSI = 0.01000 
NY =-0.00250 


NDELTA=-0.00139 

WRITE (*,1001) 

READ (*,*) WNMIN , WNMAX , [WN 
WRITE (*,1002) 

READ (#,*) XDMIN, XDMAX , IXD 
WRITE (+*,1003) 

READ (*,*) ZETA 
WRITE(* , 1100) 

READ (*,*) TL 

TL=TL*U/L 


DVR =(IZ-NRDOT) * (MASS-YVDOT) - 
& (MASS*XG-YRDOT) * (MASS*XG-NVDOT) 
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& 


AA11=((IZ-NRDOT) * YPSI-(MASS*XG-YRDOT) *NPSI) /DVR 

AA12=( (IZ-NRDOT) *YV-(MASS*XG-YRDOT) *NV) /DVR 

AA21=((NVDOT-MASS*XG) *YPSI+(MASS-YVDOT) *NPSI) /DVR 

AA22=( (MASS-YVDOT) *NV-(MASS*XG-NVDOT) *YV) /DVR 

AA13=( CIZ-NRDOT) * (YR-MASS) + 
(YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 

AA23=( (NVDOT-MASS*XG) * (YR-MASS) + 
(MASS-YVDOT) * (NR-MASS*XG) ) /DVR 

AA14=( (IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 

AA24=((NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 


BB1 =((IZ-NRDOT) * YDELTA-(MASS*XG-YRDOT) *NDELTA) /DVR 
BB2 =((MASS-YVDOT) *NDELTA-(MASS*XG-NVDOT) *YDELTA) /DVR 


ANOM=AA23 
BNOM=BB2 
CNOM=AA21 


EPS »=12D-5 
ILMAX=1500 


DO 1 I=1,IWN 
WRITE (*,2001) I,IWN 
WN=WNMIN+ (I-1) * (WNMAX-WNMIN) / CIWN-1) 
Ki=- (WN**2.DO/BNOM) - (CNOM/BNOM) 
K2=- (ANOM+2.D0*ZETA*WN) /BNOM 
DO 2 J=1,1XD 
XD=XDMIN+(J-1) * (XDMAX-XDMIN) / (IXD-1) 


CALL LINEAR(K1i,K2,XD,AA11,AA12,AA13,AA14, 
AA21,AA22,AA23,AA24,BB1i1,BB2,A,B,TL) 

CALL RGG(5,5,A,B,ALFR, ALFI,BETA,0,ZZZ,IERR) 
DO 11 IJE=1,5 

WR CIJE)=ALFR (IJE) /BETA(IJE) 

WI (IJE)=ALFI (IJE) /BETA(IJE) 
CONTINUE 
CALL DSTABL(DEOS ,WR,WI,FREQ) 


IF VC4Gi come sc 
DEOSOOG=DEOS 

XDOO =XD 

LL=0 

GOlTEr2 

DEOSNN=DEOS 

XDNN =XD 
PR=DEOSNN*DEOSOO 

IF (PR.GT.0.DO) GO TO 3 
LL=LL+1 


06 


2 


ie SCLeGr.3) STOP 1000 
IL=0 

XDO=XDOO 
XDN=XDNN 
DEOSO=DEOSO00 
DEOSN=DEOSNN 
XDL=XDO 

XDR=XDN 
DEOSL=DE0SO 
DEOSR=DEOSN 
XD=(XDL+XDR) /2.DO 


CALL LINEAR(Ki,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22, 
AA23,AA24,BBi,BB2,A,B,TL) 
CALL RGG(5,5,A,B,ALFR,ALFI,BETA,0,2ZZZ,IERR) 
DO 12 IJE=1,5 
WR(IJE)=ALFR(IJE) /BETA(IJE) 
WI (IJE)=ALFI (IJE) /BETA(IJE) 
CONTINUE 
CALL DSTABL(DEOS, WR, WI,FREQ) 


DEOSM=DEO0S 

XDM=XD 

PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 

IF (PRL.GT.0.DO) GO TO 5 
XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDL-XDM) 

IF (DIF.GT.EPS) GO TO 6 
XD=XDM 

GO TO 4 

HPO CPRReG!.0.D0) STOP 3200 
XDO=XDM 

XDN=XDR 

DEOSO=DE0SM 

DEOSN=DEOSR 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDM-XDR) 

If (DIF-GT-.EPS) C0) T0m6 
XD=XDM 

LLL=10+LL 

WRITE (LLL,*) XD,WN 


of 


3 


2 
i! 


1001 
1002 
1003 
1100 
2001 


XDOO=XDNN 
DEOSOO=DEOSNN 
CONTINUE 
CONTINUE 


FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF WN (dimensionless) ’) 
FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF XD (dimensionless) ’) 
FORMAT (’ ENTER DAMPING RATIO’) 

FORMAT (’ ENTER TIME LAG TL (dimensional) ’) 

FORMAT (215) 

END 


SUBROUTINE DSTABL(DEOS,WR,WI,OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(S) ,WI(5) 
DEOS=-1.0D+20 
DO 1 I=1,5 
IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR (I) 
I J=I 
CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 


SUBROUTINE LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23, AA24 
,BB1,BE2 Ase aL) 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION K1,K2 

DIMENSION A(5,5),B(5,5) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 

A(1,5)=0.0D0 

A(2,1)=0.0D0 

A(2,2)=0.0D0 

A(2,3)=0.0D0 

A(2,4)=0.0D0 

A(2,5)=1.0D0 
A(3,1)=AA11-(BB1*K1*TL/XD) +BB1*K1 
A(3,2)=AA12-BB1*K1*TL/XD 
A(3,3)=AA13+BB1*K2+ (BB1*K1*TL*TL/ (2.D0*XD)) 
A(3,4)=AA14+BB1*K1/XD 

A(3,5)=-1.D0+ (BB1*K1*TL*TL/ (2.D0*XD) ) 
A(4,1)=1.0D0 
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A(4,2)=1.0D0 

A(4,3)=0.0D0 

A(4,4)=0.0D0 

A(4,5)=0.0D0 
A(5,1)=AA21+BB2*K1-BB2*K1*TL/XD 

A(5, 2) =AA22-BB2*K1*TL/XD 

A(5, 3) =AA23+BB2*K2+ (BB2*K1*TL*TL/ (2.D0*XD)) 
A(5, 4) =AA24+BB2*K1/XD 
A(5,5)=(BB2*K1*TL*TL/ (2.D0*XD) ) 


B(1,1)=1.0D0 

B(1,2)=0.0D0 

B(1,3)=0.0D0 

B(1,4)=0.0D0 

B(1,5)=0.0D0 

B(2,1)=0.0D0 

B(2,2)=1.0D0 

B(2,3)=0.0D0 

B(2,4)=0.0D0 

B(2,5)=0.0D0 

—-B(3,1)=0.0D0 

B(3,2)=0.0D0 
B(3,3)=(BB1*K1*TL*TL*TL/ (6 .DO*XD) ) 
B(3,4)=0.0D0 
B(3,5)=(BB1*K1*TL*TL*TL/ (6 .D0O*XD) ) 
B(4,1)=0.0D0 

B(4,2)=0.0D0 

B(4,3)=0.0D0 

B(4,4)=1.0D0 

B(4,5)=0.0D0 

B(5,1)=0.0D0 

B(5,2)=0.0D0 

B(5,3)=1.0D0+ (BB2*K1*TL*TL*TL/ (6.D0*XD) ) 
B(5,4)=0.0D0 
B(5,5)=(BB2*K1*TL*TL*TL/ (6 .DO*XD) ) 
RETURN 

END 


©) 


QO 


PROGRAM THESIS4.FOR (Time Delay-ist Order Approx in delta ,T1) 
BIFURCATION ANALYSIS 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION K1i,K2,L,NR,NV,IZ,MASS,NDELTA,NPSI,NY, 


& 


NRDOT , NVDOT 


DIMENSION A(4,4) ,B(4,4) ,FV1(4) ,1V1(4) ,ZZZ(4,4) , ALFR(4) , 


& 


ALFI (4) ,BETA(4) ,WR(4) , W1I(4) 


OPEN (11,FILE=’BIF1.RES’) 
OPEN (12,FILE=’BIF2.RES’) 
OPEN (13,FILE=’BIF3.RES’) 


Vehicle Parameters: 


IZ =0.0 

Ie =528 

RHO =1.94 

G =32.2 

XG =0.0 

MASS =0.0088 

U =3.0 
YRDOT = 0.00000 
YVDOT =-0.00912 
YR =+0.00456 
YV =-0.01434 
YPSI = 0.01400 
ae = 0.02000 
YDELTA= 0.00278 
NRDOT =-0.00115 
NVDOT = 0.00000 
NR =-0 .00296 
NV =-0 .00460 
NPSI = 0.01000 
NY = -0.00250 


NDELTA=-0.00139 


WRITE (*,1001) 


READ (*,*) WNMIN , WNMAX, IWN 
WRITE (*,1002) 
READ (*,*) XDMIN , XDMAX , IXD 


WRITE (*,1003) 
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11 


& 


& 


& 


READ (*,*) ZETA 
WRITE(*, 1100) 

READ (*,*) TL 
TL=TL*U/L 


DVR =(IZ-NRDOT) * (MASS-YVDOT) - 
(MASS*XG-YRDOT ) * (MASS*XG-NVDOT) 


AA11=( (IZ-NRDOT) *YPSI-(MASS*XG-YRDOT) *NPSI) /DVR 
AA12=((IZ-NRDOT ) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21=( (NVDOT-MASS*XG) *YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22=( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13=((IZ-NRDOT) * (YR-MASS) + 
(YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23=((NVDOT-MASS*XG) * (YR-MASS) + 
(MASS-YVDOT) * (NR-MASS*XG) ) /DVR 
AA14=( (IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24=( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 


BB1 =((IZ-NRDOT) *YDELTA- (MASS *XG-YRDOT) *NDELTA) /DVR 
BB2 =((MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 


ANOM=AA23 
BNOM=BB2 


EPS) =1.D-5 
ILMAX=1500 


DO 1 I=1,IWN 
WRITE (*,2001) I,IWN 
WN=WNMIN+ (I-1) * (WNMAX-WNMIN) / CIWN-1) 
Ki=-WN**2.D0/BNOM 
K2=- (ANOM+2.DO*ZETA*WN) /BNOM 
DO 2 J=1,1XD 
XD=XDMIN+(J-1) * (XDMAX-XDMIN) / (IXD-1) 


CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14, 
AA21,AA22,AA23,AA24,BB1,BB2,A,B,TL) 

CALL RGG(4,4,A,B,ALFR,ALFI,BETA,0,ZZZ, IERR) 
DO 11 IJE=1,4 

WR CIJE)=ALFR (IJE) /BETA CIJE) 

WI (IJE)=ALFI (IJE) /BETA (CIJE) 
CONTINUE 
CALL DSTABL(DEOS, WR, WI ,FREQ) 


Tree cr.1) GO TO 10 


DEOSOO=DEOS 
XDOO =XD 


61 


10 


eZ 


LL=0 

GOPTURZ 

DEOSNN=DEOS 

XDNN =XD 
PR=DEOSNN*DEOSOO 

IF (PR.GT.0.D0) GO TO 3 
LL=LL+1 

IF (LL.GT.3) STOP 1000 
IL=0 

XDO=XDO0O0 

XDN=XDNN 

DEOSO=DEOSO00 
DEOSN=DEOSNN 

XDL=XDO 

XDR=XDN 

DEOSL=DEOSO 
DEOSR=DEOSN 
XD=(XDL+XDR) /2.D0 


CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22, 
AA23,AA24,BB1,BB2,A,B,TL) 
CALL RGG(4,4,A,B,ALFR,ALFI,BETA,0,ZZZ,IERR) 
DO 12 IJE=1,4 
WRCIJE)=ALFR (CIJE) /BETA(CIJE) 
WI (IJE) =ALFI (IJE) /BETA (IJE) 
CONTINUE 
CALL DSTABL(DEOS,WR,WI,FREQ) 


DEOSM=DEOS 

XDM=XD 

PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 

IF (PRL.GT.O.DO) GO TO 5 
XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDL-XDM) 

TF (DiF GT -EPS) Go T0n6 
XD=XDM 

GO TO 4 

IF (PRR.GT.0.DO) STOP 3200 
XDO=XDM 

XDN=XDR 

DEOSO=DEOSM 

DEOSN=DEOSR 
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IL=IL+1 
IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDM-XDR) 
IF (DIF.GT.EPS) GO TO 6 
XD=XDM 
LLL=10+LL 
WRITE (LLL,*) XD,WN 
XDOO=XDNN 
DEOSOO=DEOSNN 

CONTINUE 

CONTINUE 


FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF WN (dimensionless) ’) 
FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF XD (dimensionless) ’) 
FORMAT (’ ENTER DAMPING RATIO’) 

FORMAT (’? ENTER TIME LAG TL (dimensional) ’) 

FORMAT (215) 

END 


SUBROUTINE DSTABL(DEOS,WR, WI ,OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4) ,WI(4) 
DEOS=-1 .0D+20 
DO 1 I=1,4 
TE SVR CL) ei. DEOS) GO TO 1 
DEOS=WR (1) 
IJ=I 
CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 


SUBROUTINE LINEAR (K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23,AA24 
BEd Bebo. AB. Tb) 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION K1,K2 

DIMENSION A(4,4) ,B(4,4) 
A(1,1)=0.0DO 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 


A(2,1)=AA11+BB1*K1-BB1*K1*TL/XD 
A(2,2)=AA12-BB1*K1*TL/XD 
A(2,3)=AA13+BB1*K2-BB1i*K1*TL 
A(2,4)=AA14+BB1*K1i/XD 
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A(3,1)=AA21+BB2*K1-BB2*K1*TL/XD 
A(3,2)=AA22-BB2*K1*TL/XD 
A(3,3)=AA23+BB2*K2-BB2*K1*TL 
A(3,4)=AA24+BB2*K1/XD 


A(4,1)=1.0D0 
A(4,2)=1.0D0 
A(4,3)=0.0D0 
A(4,4)=0.0D0 


B(1,1)=1.0D0 
B(1,2)=0.0D0 
B(1,3)=0.0D0 
B(1,4)=0.0D0 


B(2,1)=0.0D0 
B(2,2)=1.0D0 
B(2,3)=BB1*K2*TL 
B(2,4)=0.0D0 


B(3,1)=0.0D0 
B(3,2)=0.0D0 
B(3,3)=1.0D0+ (BB2*K2*TL) 
B(3,4)=0.0D0 


B(4,1)=0.0D0 
B(4,2)=0.0D0 
B(4,3)=0.0D0 
B(4,4)=1.0D0 


RETURN 
END 


OQ 


Q 


PROGRAM THESISS.FOR (Time Delay-ist Order Approx in delta & y ie. in Ti & T2)) 
BIFURCATION ANALYSIS 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION Ki,K2,L,NR,NV,IZ,MASS,NDELTA,NPSI,NY, 


& NRDOT , NVDOT 
DIMENSION A(4,4),B(4,4) ,FV1(4) ,1V1(4) ,ZZZ(4,4) , ALFR(4), 
& ALFI (4) ,BETA(4) , WR(4) , WI (4) 


OPEN (11,FILE=’BIF1.RES’) 
OPEN (12,FILE=’BIF2.RES’ ) 
OPEN (13,FILE=’BIF3.RES’) 


Vehicle Parameters: 


IZ =0.0 

L =528 
RHO =1.94 

G =32.2 
XG =0.0 
MASS =0.0088 
U =3.0 


YRDOT = 0.00000 
YVDOT =-0.00912 
YR =+0.00456 
YV =-0.01434 


YPSI = 0.01400 
02000 


2 


YDELTA= 0.00278 
NRDOT =-0.00115 
NVDOT = 0.00000 


NR =-0 .00296 
NV =-0.00460 
NPSI = 0.01000 
NY = -0.00250 


NDELTA= -0.00139 


WRITE (*,1001) 

READ (*,*) WNMIN , WNMAX, IWN 
WRITE (+*,1002) 

READ (*,*) XDMIN , XDMAX , IXD 
WRITE (*,1003) 

READ (*,*) ZETA 
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& 


& 


11 


WRITE(* , 1100) 


READ (*,*) Thi 
TLI=TLis0/ Lc 
WRITE(*,1101) 

READ (*,*) Tee 


TL2=TL2*U/L 


DVR =(IZ-NRDOT) * (MASS-YVDOT) - 
(MASS*XG-YRDOT) * (MASS*XG-NVDOT) 


AA11=( (IZ-NRDOT) *YPSI- (MASS*XG-YRDOT) *NPSI) /DVR 
AA12=( (IZ-NRDOT) *YV- (MASS*XG-YRDOT) *NV) /DVR 
AA21=( (NVDOT-MASS*XG) * YPSI+(MASS-YVDOT) *NPSI) /DVR 
AA22=( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13=( CIZ-NRDOT) * (YR-MASS) + 
CYRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23=( (NVDOT-MASS*XG) * (YR-MASS) + 
(MASS-YVDOT) * (NR-MASS*XG) ) /DVR 
AA14=((IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24=( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 


BB1 =((IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2 =((MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA)/DVR 


ANOM=AA23 
BNOM=BB2 
CNOM=AA21 


EPS =1-D-5 
ILMAX=1500 


DO 1 I=1,IWN 
WRITE (*,2001) I,IWN 
WN=WNMIN+ (I-1) * CWNMAX-WNMIN) / C(IWN-1) 
Ki=- (WN**2.D0/BNOM) -(CNOM/BNOM) 
K2=- (ANOM+2.D0*ZETA*WN) /BNOM 
DOT 2 J=1, 
XD=XDMIN+ (J-1) * (XDMAX-XDMIN) / (IXD-1) 


CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14, 


& AA21,AA22,AA23,AA24,BB1,BB2,A,B,TL1,TL2) 


CALL RGG(4,4,A,B,ALFR,ALFI,BETA,0,ZZZ,IERR) 
DO 11 IJE=1,4 
WR(IJE)=ALFR(IJE) /BETA(IJE) 
WI(IJE)=ALFI(IJE) /BETA(IJE) 
CONTINUE 
CALL DSTABL(DEOS,WR,WI,FREQ) 
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10 


12 


iF (J.GT.1) (GO TO 10 
DEOSOO=DE0S 

XDO0O0 =XD 

LL=0 

GO TO 2 

DEOSNN=DEOS 

XDNN =XD 
PR=DEOSNN*DEOSOO 

IF (PR.GT.0.D0) GO TO 3 
LL=LL+1 

ie (eEeGr. 3) STOP 1000 
IL=0 

XDO=XDO00 

XDN=XDNN 
DEOSO=DEO0SO00 
DEOSN=DEOSNN 

XDL=XDO 

XDR=XDN 

DEOSL=DE0SO 
DEOSR=DEOSN 

XD= (XDL+XDR) /2.D0 


CALL LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14,AA21,AA22, 
AA23 ,AA24,BB1,BB2,A,B,TL1,TL2) 
CALL RGG(4,4,A,B,ALFR,ALFI,BETA,0,ZZZ, IERR) 
DO 12 IJE=1,4 
WR CIJE) =ALFR (IJE) /BETA (IJE) 
WI (IJE)=ALFI (IJE) /BETA(IJE) 
CONTINUE 
CALL DSTABL(DEOS,WR,WI,FREQ) 


DEOSM=DEO0S 

XDM=XD 

PRL=DEOSL*DEOSM 
PRR=DEOSR*DEOSM 

IF (PRL.GT.0.D0) GO TO 5S 
XDO=XDL 

XDN=XDM 

DEOSO=DEOSL 

DEOSN=DEOSM 

IL=IL+1 

IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDL-XDM) 

EDIE ]Gi. EPs) GO) 106 
XD=XDM 

GO TO 4 

IF (PRR.GT.0.DO) STOP 3200 
XDO=XDM 
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2 
it 


1001 
1002 
1003 
1100 
1101 
2001 


& 


XDN=XDR 
DEOSO=DEOSM 
DEOSN=DEOSR 
IL=IL+1 
IF (IL.GT.ILMAX) STOP 3100 
DIF=DABS (XDM-XDR) 
IF (DIF.GT.EPS) GO TO 6 
XD=XDM 
LLL=10+LL 
WRITE (LLL,*) XD,WN 
XDOO=XDNN 
DEOSOO=DEOSNN 
CONTINUE 
CONTINUE 


FORMAT (’? ENTER MIN, MAX, AND INCREMENTS OF WN (dimensionless) ’) 
FORMAT (? ENTER MIN, MAX, AND INCREMENTS OF XD (dimensionless) ’) 
FORMAT (? ENTER DAMPING RATIO’) 

FORMAT (’? ENTER TIME LAG TL1 (dimensional) ’) 

FORMAT (’ ENTER TIME LAG TL2 (dimensional) ’) 

FORMAT (215) 

END 


SUBROUTINE DSTABL(DEOS, WR, WI, OMEGA) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4) ,WI(4) 
DEOQS=-1 .0D+20 
DO 1 I=1,4 
IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR(I) 
hia 
CONTINUE 
OMEGA=WI(IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 


SUBROUTINE LINEAR(K1,K2,XD,AA11,AA12,AA13 ,AA14,AA21,AA22, AA23,AA24 
,BB1, BBQ, AB thi 2) 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION K1,kK2 

DIMENSION A(4,4) ,B(4,4) 
ACL, 1) =9,0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0D0 
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A(2,1)=AA11+BB1*K1-BB1*#K1*TL2/XD-BB1*K1*TL1i/XD 
A(2,2)=AA12-BB1*K1*TL2/XD-BB1*K1*TL1/XD 
A(2,3)=AA13+BB1*K2-BB1*K1*TL1+BB1*K1*TL1*TL2/XD 
A(2,4)=AA14+BB1*K1/XD 


A(3,1)=AA21+BB2*K1-BB2*K1*TL2/XD-BB2*K1*TL1/XD 
A(3,2)=AA22-BB2*K1*TL2/XD-BB2*K1*TL1/XD 

A(3,3) =AA23+BB2*K2-BB2*K1*TL1+BB2*K1*TL1*TL2/XD 
A(3,4) =AA24+BB2*K1/XD 


A(4,1)=1.0D0 
A(4,2)=1.0D0 
A(4,3)=0.0D0 
A(4,4)=0.0D0 


B(1,1)=1.0D0 
B(1,2)=0.0D0 
B(1,3)=0.0D0 
B(1,4)=0.0D0 


B(2,1)=0.0D0 

B(2,2)=1.0D0- (BB1*K1*TL1*TL2/XD) 
B(2,3)=BB1*K2*TL1 

B(2,4)=0.0D0 


B(3,1)=0.0DO 
B(3,2)=-BB2*K1*TL1*TL2/XD 
B(3,3)=1.0D0+(BB2*K2*TL1) 
B(3,4)=0.0D0 


B(4,1)=0.0D0 
B(4,2)=0.0D0 
B(4,3)=0.0D0 
B(4,4)=1.0D0 


RETURN 
END 
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PROGRAM HOPF.FOR 
Hopf Bifurcation Analysis 
Third Order Expansions: First Order Approximation 


IMPLICIT DOUBLE PRECISION she Oe) 

DOUBLE PRECISION K1i,K2,K3,L,NR,NV,NDELTA,NPSI,NY,IZ,MASS, 
NRDOT,NVDOT,KiP,K2P,NVVV,NVRR,NRVV,NYYY,NPPP 

DOUBLE PRECISION M11,M12,M13,M14,M21,M22,M23,M24, 
M31,M32,M33,M34,M41,M42 ,M43,M44, 
N11,N12,N13,N14,N21,N22 ,N23,N24, 
N31,N32,N33,N34,N41,N42,N43,N44, 
L21,L22,L23,124, L3a9bs2,b3s, ss, 
L41,L42 ,L43,L44 


DIMENSION A(4,4) ,T(4,4) , TINV(4,4) ,FV1(4) ,IV1(4) ,ZZZ(4, 4) 


DIMENSION WR(4) ,WI(4) , TSAVE(4,4) , TLUD(4,4) , IVLUD(4) ,SVLUD (4) 
DIMENSION ASAVE(4,4) 


OPEN (11,FILE=’BIF1.RES’ ) 
OPEN (15,FILE=’HOPF.RES’ ) 


Vehicle Parameters: 


EZ =0.0 

L =528 
RHO =1.94 

G =32.2 
XG =0.0 
MASS =0.0088 
U =1.0 


YRDOT = 0.00000 
YVDOT =-0.00912 
YR =+0.00456 
YV =-0.01434 
YVVV  =-0.15391 
YVRR =-0.05476 


YRVV = 0.04608 
YYYY = 0.46800 
YPPP = 0.00000 
YPSI = 0.01400 
er = 0.02000 


YDELTA=+0 .00278 
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NRDOT =-0.00115 
NVDOT = 0.00000 
NR =-0 .00296 
NV =-0.00460 


NVVV =-0.00336 
NVRR = 0.00759 
NRVV =-0.05160 
NYYY = 0.00000 
NPPP = 0.00000 
NPSI = 0.01000 
NY =-0 .00250 
NDELTA=-0 .00139 


WRITE (*, 1003) 

READ (*,*) ZETA 
WRITE(*,1100) 

READ (*,*) TL 
TL=TL*U/L 

WRITE (*,1006) 

READ (*,*) DO 


DVR =(IZ-NRDOT) * (MASS-YVDOT) - 
& (MASS*XG-YRDOT) + (MASS*XG-NVDOT ) 


AA11=((IZ-NRDOT) *YPSI-(MASS*XG-YRDOT) *NPSI) /DVR 
AA12=( (IZ-NRDOT) *YV-(MASS*XG-YRDOT) *«NV) /DVR 
AA21=( (NVDOT-MASS*XG) * YPSI+ (MASS-YVDOT) *NPSI) /DVR 
AA22=( (MASS-YVDOT) *NV- (MASS*XG-NVDOT) *YV) /DVR 
AA13=( (IZ-NRDOT) * (YR-MASS) + 


& (YRDOT-MASS*XG) * (NR-MASS*XG) ) /DVR 
AA23=( (NVDOT-MASS*XG) * (YR-MASS) + 
& (MASS-YVDOT) * (NR-MASS*XG) ) /DVR 


AA14=((IZ-NRDOT) *YY+ (YRDOT-MASS*XG) *NY) /DVR 
AA24=( (NVDOT-MASS*XG) *YY+ (MASS-YVDOT) *NY) /DVR 


BB1 =((IZ-NRDOT) *YDELTA- (MASS*XG-YRDOT) *NDELTA) /DVR 
BB2 =((MASS-YVDOT) *NDELTA- (MASS*XG-NVDOT) *YDELTA) /DVR 


ANOM=AA23 
BNOM=BB2 
CNOM=AA21 


EPS =1.D=5 
ILMAX=1500 


IWN=1000 
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DO 1 II=1,IWN 


WRITE (*,2001) II 

READ (11,*) XD,WN 

K1=- (WN**2.D0/BNOM) - (CNOM/BNOM) 
K2=- (ANOM+2 .D0O*ZETA*WN) /BNOM 


K3=K2 


Start Hopf Bifurcation Analysis 


Evaluate Nonlinear Rudder Expansion Coefficients 


K2=0 . ODO 


KiP=Ki- 
K2P=K2- 


DPPV=- 

z 
DPVV=- 

- 
DPPR=- 
DPRR=- 
DPPY=- 
DPYY=- 

_ 
DVVR=- 
DVRR=- 
DVVY=- 
DVYY=- 

5 
DRRY=- 
DRYY=- 
DPVR=- 
DPVY=- 
DPRY=- 
DVRY=- 
DPPP=- 

5 
DVVV=- 

: 
DRRR=- 
DYYY=- 


K1*TL*U/XD 
K1*TL/XD 


(1.D0/ (3.D0*DO**2)) *3.D0*K1P*K1P*K2P 
0.5D0*K1*TL/XD + 3.D0*K1* (TL*U) **3/ (3 .DO*XD**3) 
(1.D0/(3.D0*D0**2) ) *3.DO+#K1P+*K2P+*K2P 
3.DO*U*TL*TL*TL*K1/ (3.D0*XD**3) 
(1.D0/(3.D0*DO**2) )*3.D0*K1P*K1P*K3 
(1.D0/(3.DO*DO**2) ) *3..DO*K1P*K3*K3 

(1.D0/ (3 .DO*D0O**2) ) *3.DO*K1P*K1P*K1/XD 
3.D0O*TL*TL*U*U*K1/ (3 .DO*XD**3) 
(1.D0/(3.DO*D0**2) ) *3.DO*K1P*K1*K1/ (XD**2) 
3.DO*#TL*U*K1/ (3 .D0*XD**3) 
(1.D0/(3.D0*D0**2) ) *3.D0*K2P*K2P*K3 
(1.D0/(3.D0*D0**2)) *3.D0*K2P*K3*K3 
(1.D0/(3.D0*D0**2)) *3.D0*K1*K2P*K2P/XD 
3.D0*K1*TL*TL/ (3.D0*XD**3) 

(1.D0/ (3 .D0*DO**2) ) *3.DO*K1*K1*K2P/ (XD**2) 
3.D0*TL*K1/(3.D0*XD**3) 

(1.D0/(3.D0*D0**2) )*3.D0*K1*K3*K3/XD 
(1.D0/(3.D0*#DO*#*2) ) *3.D0#*#K1*K1*K3/ (XD**2) 
(1.D0/(3.D0*DO**2) )*6.D0*K1P*K2P*K3 
(1.D0/(3.D0*DO**2) ) *6.DO*K1P*K1*K2P/XD 

6 .DO*TL*TL*U*K1/ (3. DO*XD**3) 
(1.D0/(3.D0*DO**2) ) *6.D0*K1P*K1*K3/XD 
(1.D0/(3.D0*D0**2) ) *6 .DO*K1*K2P*K3/XD 

(1.D0/ (3 .D0*D0**2) ) *1.D0*K1P*K1P*K1P 
Ki*TL*U/(6.D0*XD) + (K1* (TL*U) **3)/(3.D0*XD**3) 
(1.D0/(3.DO*DO**2) ) *1.D0*K2P*K2P*K2P 

K1* (TL**3) /(3.D0*XD**3) 
(1.D0/(3.D0*DO**2) ) *1.D0*K3*K3*K3 

(1.DO0/ (3.D0*D0**2) ) * (K1/XD) **3-K1/ (3.D0*XD**3) 
K1/ (3.D0*XD**3) 


Evaluate Transformation Matrix of Eigenvectors 


@Z 


OQ 


OQ 


i 
it 


30 


14 


60 


15 


70 


13 


16 
ee 


CALL LINEAR(K1,K3,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23, 
AA24,BB1,BB2,A,TL) 


DO 11 I1=1,4 
DO 12 J=1,4 
ASAVE(I,J)=A(I, J) 
CONTINUE 
CONTINUE 
CALL RG(4,4,A,WR,WI,1,ZZZ,IV1,FV1,IERR) 
CALL DSOMEG(IEV,WR,WI,OMEGA , CHECK) 
OMEGAO=OMEGA 
DO 50 1=1,4 
Tiel, 1) = 2227. EeLEV) 
mG, 2) =-7ZZ. Clee EV +1 ) 
CONTINUE 
IF (CIEV-EQ ive COeTO 13 
IF (IEV.EQ.2) GO TO 14 
TF CIEV EQs3 eco TO 15 
STOP 3004 
DO 60 I=1,4 
TGS) —=22201. 
TCI ,4)=ZZZ(1,4) 
CONTINUE 
GO TO 17 
DO 70 I1=1,4 
Tis) —Zerl, 1) 
T(1,4)=ZZZ(1,2) 
CONTINUE 
GO TO 17 
DO 16 I=1,4 
TGs) —ZZZ2( 1,3) 
TCI ,4)=Z2ZZ(1, 4) 
CONTINUE 
CONTINUE 


Normalization of the Critical Eigenvector 
CALL NORMAL (T) 

Definition of Mij 

Mii=T(1,1) 

M21=T(2,1) 

M31=T(3, 1) 

M41=T (4,1) 


M22=T (2,2) 
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M32=T (3,2) 
M42=T (4,2) 
M13=T(1,3) 
M23=T (2,3) 
M33=T(3,3) 
M43=T (4,3) 
M14=T(1,4) 
M24=T (2,4) 
M34=T (3,4) 
M44=T (4,4) 


Definition of Lij 


L21= DPPV*M11*M11*M21 + DPVV*M11*M21*M21 + DPPR*M11*M11*M31 
+ DPRR*M11*M31*M31 + DPPY*M11*M11*M41 + DPYY*M11*M41*M41 
+ DVVR*M21*M21*M31 + DVRR*M21*M31*M31 + DVVY*M21*M21+*M41 
+ DVYY*M21*M41*M41 + DRRY*M31*M31*M41 + DRYY*M31*M41*M41 
+ DPVR*M11*M21*M31 + DPVY*M11*M21*M41 + DPRY*M11*M41*M31 
+ DVRY*M21*M41*M31 + DPPP*M11*M11*M11 + DVVV*M21*M21*M21 
+ DRRR*M31*M31*M31 + DYYY*M41*M41*M41 

PPV = M11*M114*#M22 + 2.0*#M11*M12*M21 

PVV M1i2*M21*M21 + 2.0*M11*M21*M22 

PPR = Mi1*M11*M32 + 2.0*M11*M12*M31 

PRR = M12*M314*M31 + 2.0*M11*M31*M32 

PPY = M11*M11*M42 + 2.0*M11*M12+*M41 

PYY = M41*M41*M1i2 + 2.0*M11*M41+*M42 

VVR = M21*M21*M32 + 2.0*M31*M21*M22 

VRR = M22*M31*M31 + 2.0*M31*M32*M21 

VVY = M21*M21*M42 + 2.0*M41*M21*M22 

VYY = M22*M41*M41 + 2.0*M41*M42*M21 

RRY = M31*M31*M42 + 2.0*M41*M31*M32 

RYY = M32*M41*M41 + 2.0*M41*M42*M31 

PVR = M11*M21*M32+M31* (M11*M22+M12*M21) 

PVY = M11*M21*M42+M41* (M11*M22+M12*M21) 

PRY = M11*M41*M32+M31* (M11*M42+M12+*M41) 

VRY = M21*M41*M32+M31* (M21*M42+M22+*M41) 

PPP = 3.0*#M11*«M11*M12 

VV = 3.0*M21*M21*M22 

RRR = 3.0*M31*M31*M32 

YYY = 3.0*M41*M41*M42 


L22=DPPV*PPV+DPVV*PVV+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
+DVVR*VVR+DVRR*VRR+DVVY*VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 
+DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVVV*VVV 
+DRRR*RRR+DYYY*YYY 


PPV = M12*M12*M21 + 2.0*M11*M124*M22 
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RP BR RB BR kB & 


M11*M22*M22 
M12*M12*M31 
M11*M32*M32 
M12*M12*M41 
M11*M42*M42 
M22*M22*M31 
M21*M32*M32 


= M22*M22*M41 
= M21*M42+*M42 
= M32*M32*M41 


M31*M42*M42 


= M12*M22*M31 


M12*M22*M41 
M12*M42*M31 


= M22*M42+*M31 
3.0*M11*M12*M12 
3 .0*M21*M22*M22 
3 .0*M31*M32*M32 
3 .0*M41*M42*M42 


t+ ttt tte tt t+ t+ + + t+ 


+ 


.0*M124*M21*M22 
.O*M11*M12*M32 
.0*M12*M31*M32 
.0O*M11*M12*M42 
.0*M12*M41*M42 
.0*M21*M22*M32 
. O#¥M22*M314*M32 
. O*M214*M22*M42 
. O*M22*M41*M42 
. O*M31*M32*M42 
.0*M32*M41*M42 
M32* (M11*M22+M12*M21) 
M42* (M11*M22+M12*M21) 
M32* (M11*M42+M12*M41) 
M32* (M21*M42+M22*M41 ) 


NONONNNNNNN ND WN 


L23=DPPV*PPV+DPVV*PVV+DPPR*PPR+DPRR*PRR+DPPY*PPY+DPYY*PYY 
+DVVR*VVR+DVRR*VRR+DVVY*VVY+DVYY*VYY+DRRY*RRY+DRYY*RYY 
+DPVR*PVR+DPVY*PVY+DPRY*PRY+DVRY*VRY+DPPP*PPP+DVVV*VVV 


+DRRR*RRR+DYYY*YYY 

L24= DPPV*M12*M12*M22 + DPVV*M12*M22*M22 
+ DPRR*M12*M32*M32 + DPPY*M12*M12*M42 
+ DVVR*M22*M22*M32 + DVRR*M22*M32*M32 
+ DVYY*M22+M42*M42 + DRRY*M32*M32*M42 
+ DPVR*M12*M22*M32 + DPVY*M12*M22*M42 
+ DVRY*M22*M42*M32 + DPPP*M12*M12*M12 
+ DRRR*M32*M32*M32 + DYYY*M42*M42*M42 

Bote e2 1 

L32=L22 

L33=L23 

L34=L24 


L21=L21*BB1*U*U 
L22=L22*BB1*U*U 
L23=L23*BB1*U*U 
L24=L24*BB1*U*U 
L31=L31*BB2*U*U 
L32=L32*BB2*U*U 
L33=L33*BB2*U*U 
L34=L34*BB2*U*U 
L41=-0. 5*M11*M11* (M21+U*M11/3.0) 
L42=-M11* (M12*M21+0.5*M11*M22+0 .5*U*M12*M11) 
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++ + + + + 


DPPR*M12*M12*M32 
DPYY*M12*M42*M42 
DVVY*M22*M22*M42 
DRYY*M32*M42*M42 
DPRY*M12*M42*M32 
DVVV*M22*M22*M22 


QO 


L43=-M12* (M11*M22+0.5*M12*M21+0.5*U*M11*M12) 
L44=-0.5*M12*M12* (M22+U*M12/3 .0) 


C11=(1/6.0) * YVVV* (M21*M21*#M21)+0.5*YVRRx* (M31*M31*M21)+ 
0.5*YRVV*M21*M21*M31+(1/6.0) *YPPP*M11*M11*M11+(1/6.0) *YYYY*M41*M41*M41 


C12=(1/6.0) * YVVV* (3*M21*M21*M22)+0.5*YVRR* (M31*M31* 
M22+24*M31*M32*M21)+0.5*YRVV* (M21*M21*M32+2*M31*M21*M22)+(1/6.0) *YPPP*3*M11*M1i1*M12 
+(1/6.0) *YYYY*3*M41+*M41*M42 


C13=(1/6.0) *YVVV* (3*M21*M22*M22)+0.5*YVRR* (M21*M32* 
M32+2*M31*M32*M22)+0.5*YRVV* (M22*M22*M31+2*M32*M21*M22)+(1/6.0) *YPPP*3*M11*M12*M12 
+(1/6.0)*YYYY*3*M414*M42*M42 


C14=(1/6.0) *YVVV* (M22*M22*M22)+0.5*YVRR* (M32*M32*M22)+ 
0. 5* YRVV*M22*M22*M32+ (1/6. 0) *YPPP*M12*M12*M12+(1/6.0) *YYYY*M42*M42«M42 


C21=(1/6.0) *NVVV* (M21*M21*M21)+0.5* (NVRR*M31*M31*M21)+ 
0. 5*NRVV*M21*M21*M31+(1/6.0) *NPPP*Mi1*M11*M11+(1/6.0) *NYYY*M41*M41*M41 


C22=(1/6.0) *NVVV* (34M21*M21 *M22) +0. 5*NVRR* (M31*M31* 
M22+24*M31*M32*M21)+0.5*NRVV* (M21*M21*M32+2*M31*M21*M22)+(1/6.0) *NPPP*3*M11*M11*M12 
+(1/6.0) *NYYY*3*M41*M41*M42 


C23=(1/6.0) *NVVV* (3*M21*M22+*M22) +0 .5*NVRR* (M21*M32* 
M32+2*M31*M32*M22)+0. 5*NRVV* (M22*M22*M3 1+2*M32*M21*M22) +(1/6.0) *NPPP*3*M11*M12*M12 
+(1/6.0) *NYYY*3*M41+*M424*M42 


C24=(1/6.0) *NVVV* (M22*M22*M22)+0.5*NVRR*M32*M32*M22+ 
0. 5*NRVV*M22*M22*M32+ (1/6.0)*NPPP*M12*M12*M12+(1/6.0) *NYYY*M42*M42+*M42 


D11=(C11* (IZ-NRDOT) +C21* (YRDOT-M*XG) ) /DVR 
D12=(C12* (IZ~NRDOT) +C22* (YRDOT-M*XG) ) /DVR 
D13=(C13* (IZ~NRDOT) +C23* (YRDOT-M*XG) ) /DVR 
D14=(C14* (IZ-NRDOT) +C24* (YRDOT-M*XG) ) /DVR 


D21=(C11* (NVDOT-M*XG) +C21* (M-YVDOT) ) /DVR 
D22=(C12* (NVDOT-M*XG) +C22* (M-YVDOT) ) /DVR 
D23=(C13* (NVDOT~M*XG) +C23* (M-YVDOT) ) /DVR 
D24=(C14* (NVDOT-M*XG) +C24* (M-YVDOT) ) /DVR 


Invert Transformation Matrix 


DO 20 [=1,4 
i 


DO 30 J=1,4 
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QQ 


30 
20 


40 


90 
80 


TINV(I,J)=0.0 
TSAVE(I, J)=T(I, J) 
CONTINUE 
CONTINUE 
CALL DLUD(4,4,TSAVE,4,TLUD, IVLUD) 
DOee20n—14 
IF (IVLUD(I).EQ.0) STOP 3003 
CONTINUE 
CALL DILU(4,4,TLUD,IVLUD ,SVLUD) 
DO 80 I=1,4 
DO 90 J=1,4 
TINV(I,J)=TLUD(I, J) 
CONTINUE 
CONTINUE 


Check Inv(T) *A*T 


C27) 


IMULT=2 
IF (IMULT.EQ.1) CALL MULT(TINV, ASAVE,T) 
PeeeehuULT .EQ.0) STOP 


Definition of Nij 


N11=TINV(1,1) 
N21=TINV(2,1) 
N31=TINV(3,1) 
N41=TINV(4,1) 
N12=TINV (1,2) 
N22=TINV(2,2) 
N32=TINV (3,2) 
N42=TINV (4,2) 
N13=TINV(1,3) 
N23=TINV(2,3) 
N33=TINV(3,3) 
N43=TINV (4,3) 
N14=TINV(1 ,4) 
N24=TINV(2,4) 
N34=TINV(3,4) 
N44=TINV(4,4) 


R11=N12*L21+N13*L31+N14*L41+N12*D11+N13*D21 
R12=N12*L22+N13*L32+N14*L42+N12*D12+N13*D22 
R13=N12*L23+N13*L33+N14*L43+N12*D13+N13*D23 
R14=N12*L24+N13*L34+N14*L44+N12*D14+N13*D24 
R21=N22*L21+N23*L31+N24*L414+N22*D114+N23*D21 
R22=N22*L22+N23*L32+N24*L42+N22*D12+N23*D22 
R23=N22*L23+N23*L33+N24*L43+N22*D13+N23*D23 
R24=N22*L24+N23*L34+N24*L444N22*D14+N23*D24 
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Q 


Q 


Evaluate Alpha’ and Omega’ 


DINC=0.001 
XDR =XD+DINC 
XDL =XD-DINC 
XD =XDR 


CALL LINEAR(K1,K3,XD,AA11,AA12,AA13,AA14,AA21 ,AA22,AA23, 
& AA24,BB1,BB2,A,TL) 


CALL RG(4,4,A,WR,WI,0,ZZZ,1V1,FV1, IERR) 
CALL DSTABL(DEOS ,WR,WI, FREQ) 
ALPHR=DEOS 

OMEGR=FREQ 


XD=XDL 


CALL LINEAR(K1,K3,XD,AA11,AA12,AA13,AA14,AA21,AA22,AA23, 
& AA24,BB1,BB2,A,TL) 


CALL RG(4,4,A,WR,WI,0,ZZZ,IV1,FV1, IERR) 
CALL DSTABL(DEOS, WR ,WI ,FREQ) 
ALPHL=DEOS 

OMEGL=FREQ 


DALPHA=(ALPHR-ALPHL) / (XDR-XDL) 
DOMEGA=(OMEGR-OMEGL) / (XDR-XDL) 


Evaluation of Hopf Bifurcation Coefficients 


COEF 1=3 .0*R11+R13+R22+3.0*R24 

COEF2=3 .0*R21+R23-R12-3.0*R14 

AMU2 =-COEF1/(8.0*DALPHA) 

BETA2=0 .25*COEF 1 

TAU2 =- (COEF2-~DOMEGA*COEF1/DALPHA) / (8.0*OMEGAO) 
PER =2.0*3.1415927/0MEGAO 

PER =PER*U/L 


WRITE (15,2002) XD,WN,COEF1,DALPHA, PER 


1 CONTINUE 


1001 FORMAT (’? ENTER MIN, MAX, AND INCREMENTS OF WN (dimensionless) ’) 
1002 FORMAT (’ ENTER MIN, MAX, AND INCREMENTS OF XD (dimensionless) ’) 
1003 FORMAT (’? ENTER DAMPING RATIO’) 
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1006 FORMAT (’ ENTER DSAT ’) 

1100 FORMAT (’ ENTER TIME LAG TL (dimensional) ’) 
2001 FORMAT (215) 

2002 FORMAT (5E15.5) 


SS Ye Se ey sc ey ey ey ey ee ee ey ey ey ey ey ee ee ee ee ee ee ee ee ee ee ee ey ee ee 
—_—— a Se a ee ee ee ee ee ee ee ey ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee ee eS 


SUBROUTINE DSTABL(DEOS,WR, WI, OMEGA) 
EVALUATES THE EIGENVALUE WITH THE MAXIMUM REAL PART 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4) ,WI(4) 
DEOS=-1.0D+20 
DO 1 I=1,4 
IF (WR(I).LT.DEOS) GO TO 1 
DEOS=WR (I) 
I J=I 
1 CONTINUE 
OMEGA=WI (IJ) 
OMEGA=DABS (OMEGA) 
RETURN 
END 


SUBROUTINE LINEAR(K1,K2,XD,AA11,AA12,AA13,AA14, AA21, 
& AA22,AA23,AA24,BB1,BB2,A,TL) 


FORMS THE LINEARIZED MATRIX A (time delay ist order approximation 


IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DOUBLE PRECISION K1i,kK2 
DIMENSION A(4,4) 

A(1,1)=0.0D0 

A(1,2)=0.0D0 

A(1,3)=1.0D0 

A(1,4)=0.0DO0 
A(2,1)=AA11+BB1*K1-BB1*K1i*TL/XD 
A(2,2)=AA12-BB1*K1*TL/XD 
A(2,3)=AA13+BB1*K2 
A(2,4)=AA14+BB1*K1/XD 
A(3,1)=AA21+BB2*K1-BB2*K1*TL/XD 
A(3,2)=AA22-BB2*K1*TL/XD 
A(3,3)=AA23+BB2*K2 
A(3,4)=AA24+BB2*K1/XD 
A(4,1)=1.0D0 

A(4,2)=1.0D0 


fe 


A(4,3)=0.0D0 
A(4,4)=0.0D0 
RETURN 

END 


SUBROUTINE DSOMEG(IJK,WR,WI,OMEGA, CHECK) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION WR(4) ,WI(4) 
CHECK=-1.0D+25 
DO 1 I=1,4 
IF (WR(I).LT.CHECK) GO TO 1 
CHECK=WR(I) 
ig=i 
1 CONTINUE 
OMEGA=DABS (WI (IJ)) 
IF (WI(IJ).GT.0.DO) IJK=IJ 
IF (WI(IJ).LT.0.DO) IJK=IJ-1 
RETURN 


SUBROUTINE NORMAL(T) 
IMPLICIT DOUBLE PRECISION (A-H,0-Z) 
DIMENSION T(4,4) , TNOR(4, 4) 
CFAC=T (1,1) **2+T(1,2)**2 
IF (DABS(CFAC) .LE.(1.D-10)) STOP 4001 
TNOR(1,1)=1.D0 
TNOR(2,1)=(T(1,1)*T(2,1)+T(2,2)*T(1,2))/CFAC 
TNOR(3,1)=(T(1,1) *T(3,1)+T(3,2)*T(1,2))/CFAC 
TNOR(4,1)=(T(1,1)*T(4,1)+T(4,2) *T(1,2))/CFAC 
TNOR(1,2)=0.D0 
TNOR(2,2)=(T(2,2) *T(1,1)-T(2,1)*T(1,2))/CFAC 
TNOR (3, 2)=(T(3,2)*T(1,1)-T(3,1) *T(1,2))/CFAC 
TNOR (4,2) =(T(4, 2) *T(1,1)-T(4,1)*T(1,2))/CFAC 
DO 1 I[=1,4 
DO 2 J=1,2 
T(I,J)=TNOR(I, J) 

2 CONTINUE 

1 CONTINUE 
RETURN 
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SUBROUTINE MULT(TINV,A,T) 
IMPLICIT DOUBLE PRECISION (A-H,0O-Z) 
DIMENSION TINV(4,4) ,A(4,4) ,T(4,4) ,A1(4,4) ,A2(4,4) 
DO 1 I=1,4 
DO 2 J=1,4 
A1(1I,J)=0.DO 
A2(I,J)=0.D0 
CONTINUE 
1 CONTINUE 
DO 3 I=1,4 
Does J=1,4 
DO 5 k=1,4 
A1(1I,J)=A(I,K)*T(K,J)+A1 (1, J) 
5 CONTINUE 
CONTINUE 
3 CONTINUE 
DO 6 I=1,4 
DO 7 J=1,4 
DO 8 K=1,4 
A2(1,J)=TINV(I,K)*A1(K,J)+A2(1I, J) 
8 CONTINUE 
CONTINUE 
6 CONTINUE 
DO’ 11 f=1,4 
WRITE (*,101) (ACI,J),J=1,4) 
11 CONTINUE 
DO 12 I=1,4 
WRITE (*,101) (T(I,J) ,J=1,4) 
12 CONTINUE 
DO 10 I=1,4 
WRITE (*,101) (A2(I,J),J=1,4) 
10 CONTINUE 
WRITE (*,101) A2(1,1) 
RETURN 
101 FORMAT (4E15.5) 
END 


if 


~] 


8] 





N 
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